Merge branch 'master' of https://github.com/Z3Prover/z3

README-CMake.md

@@ -1,4 +1,4 @@
-# Z3's CMake build system
+    # Z3's CMake build system
 
 [CMake](https://cmake.org/) is a "meta build system" that reads a description
 of the project written in the ``CMakeLists.txt`` files and emits a build

src/ast/rewriter/seq_rewriter.cpp

@@ -1506,6 +1506,7 @@ br_status seq_rewriter::mk_re_opt(expr* a, expr_ref& result) {
 }
 
 br_status seq_rewriter::mk_eq_core(expr * l, expr * r, expr_ref & result) {
+    TRACE("seq", tout << mk_pp(l, m()) << " = " << mk_pp(r, m()) << "\n";);
     expr_ref_vector lhs(m()), rhs(m()), res(m());
     bool changed = false;
     if (!reduce_eq(l, r, lhs, rhs, changed)) {

src/ast/seq_decl_plugin.cpp

@@ -876,6 +876,36 @@ bool seq_decl_plugin::is_value(app* e) const {
     }
 }
 
+bool seq_decl_plugin::are_equal(app* a, app* b) const {
+    if (a == b) return true;
+    // handle concatenations
+    return false;
+}
+
+bool seq_decl_plugin::are_distinct(app* a, app* b) const {
+    if (a == b) {
+        return false;
+    }
+    if (is_app_of(a, m_family_id, OP_STRING_CONST) &&
+        is_app_of(b, m_family_id, OP_STRING_CONST)) {
+        return true;
+    }
+    if (is_app_of(a, m_family_id, OP_SEQ_UNIT) && 
+        is_app_of(b, m_family_id, OP_SEQ_UNIT)) {
+        return true;
+    }
+    if (is_app_of(a, m_family_id, OP_SEQ_EMPTY) && 
+        is_app_of(b, m_family_id, OP_SEQ_UNIT)) {
+        return true;
+    }
+    if (is_app_of(b, m_family_id, OP_SEQ_EMPTY) && 
+        is_app_of(a, m_family_id, OP_SEQ_UNIT)) {
+        return true;
+    }    
+    return false;
+}
+
+
 expr* seq_decl_plugin::get_some_value(sort* s) {
     seq_util util(*m_manager);
     if (util.is_seq(s)) {

src/ast/seq_decl_plugin.h

@@ -182,7 +182,11 @@ public:
 
     virtual bool is_value(app * e) const;
 
-    virtual bool is_unique_value(app * e) const { return is_value(e); }
+    virtual bool is_unique_value(app * e) const { return false; }
+
+    virtual bool are_equal(app* a, app* b) const;
+
+    virtual bool are_distinct(app* a, app* b) const;
 
     virtual expr * get_some_value(sort * s);
 

src/shell/lp_frontend.cpp

@@ -80,7 +80,6 @@ void run_solver(lp_params & params, char const * mps_file_name) {
     solver->settings().set_message_ostream(&std::cout);
     solver->settings().report_frequency = params.rep_freq();
     solver->settings().print_statistics = params.print_stats();
-    solver->settings().presolve_with_double_solver_for_lar = params.presolve_with_dbl();
     solver->find_maximal_solution();
 
     *(solver->settings().get_message_ostream()) << "status is " << lp_status_to_string(solver->get_status()) << std::endl;

src/smt/smt_setup.cpp

@@ -471,6 +471,10 @@ namespace smt {
         setup_r_arith();
     }
 
+    void setup::setup_QF_LIRA(static_features const& st) {
+        setup_mi_arith();
+    }
+
     void setup::setup_QF_LIA() {
         TRACE("setup", tout << "setup_QF_LIA(st)\n";);
         m_params.m_relevancy_lvl       = 0;
@@ -720,10 +724,9 @@ namespace smt {
     }
 
     void setup::setup_r_arith() {
-        m_context.register_plugin(alloc(smt::theory_mi_arith, m_manager, m_params));
+        // to disable theory lra
-        
+        // m_context.register_plugin(alloc(smt::theory_mi_arith, m_manager, m_params));        
-        // Disabled in initial commit of LRA additions
+        m_context.register_plugin(alloc(smt::theory_lra, m_manager, m_params));
-        // m_context.register_plugin(alloc(smt::theory_lra, m_manager, m_params));
     }
 
     void setup::setup_mi_arith() {
@@ -937,7 +940,9 @@ namespace smt {
         }
         
         if (st.num_theories() == 1 && is_arith(st)) {
-            if (st.m_has_real)
+            if ((st.m_has_int && st.m_has_real) || (st.m_num_non_linear != 0)) 
+                setup_QF_LIRA(st);
+            else if (st.m_has_real)
                 setup_QF_LRA(st);
             else
                 setup_QF_LIA(st);

src/smt/smt_setup.h

@@ -65,6 +65,7 @@ namespace smt {
         void setup_QF_LRA();
         void setup_QF_LRA(static_features const & st);
         void setup_QF_LIA();
+        void setup_QF_LIRA(static_features const& st);
         void setup_QF_LIA(static_features const & st);
         void setup_QF_UFLIA();
         void setup_QF_UFLIA(static_features & st);

src/smt/tactic/smt_tactic.cpp

@@ -21,6 +21,7 @@ Notes:
 #include"smt_kernel.h"
 #include"smt_params.h"
 #include"smt_params_helper.hpp"
+#include"lp_params.hpp"
 #include"rewriter_types.h"
 #include"filter_model_converter.h"
 #include"ast_util.h"
@@ -64,6 +65,10 @@ public:
         return m_params;
     }
 
+    params_ref & params() {
+        return m_params_ref;
+    }
+
     void updt_params_core(params_ref const & p) {
         m_candidate_models     = p.get_bool("candidate_models", false);
         m_fail_if_inconclusive = p.get_bool("fail_if_inconclusive", true);
@@ -73,6 +78,7 @@ public:
         TRACE("smt_tactic", tout << "updt_params: " << p << "\n";);
         updt_params_core(p);
         fparams().updt_params(p);
+        m_params_ref.copy(p);
         m_logic = p.get_sym(symbol("logic"), m_logic);
         if (m_logic != symbol::null && m_ctx) {
             m_ctx->set_logic(m_logic);
@@ -84,6 +90,7 @@ public:
         r.insert("candidate_models", CPK_BOOL, "(default: false) create candidate models even when quantifier or theory reasoning is incomplete.");
         r.insert("fail_if_inconclusive", CPK_BOOL, "(default: true) fail if found unsat (sat) for under (over) approximated goal.");
         smt_params_helper::collect_param_descrs(r);
+        lp_params::collect_param_descrs(r);
     }
 
 
@@ -112,10 +119,12 @@ public:
     struct scoped_init_ctx {
         smt_tactic & m_owner;
         smt_params   m_params; // smt-setup overwrites parameters depending on the current assertions.
+        params_ref   m_params_ref;
 
         scoped_init_ctx(smt_tactic & o, ast_manager & m):m_owner(o) {
             m_params = o.fparams();
-            smt::kernel * new_ctx = alloc(smt::kernel, m, m_params);
+            m_params_ref = o.params();
+            smt::kernel * new_ctx = alloc(smt::kernel, m, m_params, m_params_ref);
             TRACE("smt_tactic", tout << "logic: " << o.m_logic << "\n";);
             new_ctx->set_logic(o.m_logic);
             if (o.m_callback) {

src/smt/theory_lra.cpp

@@ -143,14 +143,11 @@ namespace smt {
         theory_lra&          th;
         ast_manager&         m;
         theory_arith_params& m_arith_params;
-        lp_params            m_lp_params;    // seeded from global parameters.
         arith_util           a;
 
         arith_eq_adapter     m_arith_eq_adapter;
 
         vector<rational>    m_columns;
-        int m_print_counter = 0;
-
         // temporary values kept during internalization
         struct internalize_state {
             expr_ref_vector     m_terms;                     
@@ -282,14 +279,15 @@ namespace smt {
         expr*  get_owner(theory_var v) const { return get_enode(v)->get_owner(); }        
 
         void init_solver() {
+            if (m_solver) return;
+            lp_params lp(ctx().get_params());
             m_solver = alloc(lean::lar_solver); 
             m_theory_var2var_index.reset();
             m_solver->settings().set_resource_limit(m_resource_limit);
-            m_solver->settings().simplex_strategy() = static_cast<lean::simplex_strategy_enum>(m_lp_params.simplex_strategy());
+            m_solver->settings().simplex_strategy() = static_cast<lean::simplex_strategy_enum>(lp.simplex_strategy());
-            m_solver->settings().presolve_with_double_solver_for_lar = m_lp_params.presolve_with_dbl();
             reset_variable_values();
             m_solver->settings().bound_propagation() = BP_NONE != propagation_mode();
-            m_solver->set_propagate_bounds_on_pivoted_rows_mode(m_lp_params.bprop_on_pivoted_rows());
+            m_solver->set_propagate_bounds_on_pivoted_rows_mode(lp.bprop_on_pivoted_rows());
             //m_solver->settings().set_ostream(0);
         }
 
@@ -646,9 +644,11 @@ namespace smt {
         
 
     public:
-        imp(theory_lra& th, ast_manager& m, theory_arith_params& p): 
+        imp(theory_lra& th, ast_manager& m, theory_arith_params& ap): 
-            th(th), m(m), m_arith_params(p), a(m), 
+            th(th), m(m), 
-            m_arith_eq_adapter(th, p, a),
+            m_arith_params(ap), 
+            a(m), 
+            m_arith_eq_adapter(th, ap, a),            
             m_internalize_head(0),
             m_delay_constraints(false), 
             m_delayed_terms(m),
@@ -657,8 +657,8 @@ namespace smt {
             m_assume_eq_head(0),
             m_num_conflicts(0),
             m_model_eqs(DEFAULT_HASHTABLE_INITIAL_CAPACITY, var_value_hash(*this), var_value_eq(*this)),
+            m_solver(0),
             m_resource_limit(*this) {
-            init_solver();
         }
         
         ~imp() {
@@ -666,6 +666,10 @@ namespace smt {
             std::for_each(m_internalize_states.begin(), m_internalize_states.end(), delete_proc<internalize_state>());
         }
 
+        void init(context* ctx) {
+            init_solver();
+        }
+        
         bool internalize_atom(app * atom, bool gate_ctx) {
             if (m_delay_constraints) {
                 return internalize_atom_lazy(atom, gate_ctx);
@@ -1978,7 +1982,6 @@ namespace smt {
         typedef pair_hash<obj_hash<rational>, bool_hash> value_sort_pair_hash;
         typedef map<value_sort_pair, theory_var, value_sort_pair_hash, default_eq<value_sort_pair> > value2var;
         value2var               m_fixed_var_table;
-        const lean::constraint_index null_index = static_cast<lean::constraint_index>(-1);
 
         void propagate_eqs(lean::var_index vi, lean::constraint_index ci, lean::lconstraint_kind k, lp::bound& b) {
             if (propagate_eqs()) {
@@ -2008,7 +2011,7 @@ namespace smt {
 
         bool proofs_enabled() const { return m.proofs_enabled(); }
 
-        bool use_tableau() const { return m_lp_params.simplex_strategy() < 2; }
+        bool use_tableau() const { return lp_params(ctx().get_params()).simplex_strategy() < 2; }
 
         void set_upper_bound(lean::var_index vi, lean::constraint_index ci, rational const& v) { set_bound(vi, ci, v, false);  }
 
@@ -2022,10 +2025,10 @@ namespace smt {
             lean::var_index ti = m_solver->adjust_term_index(vi);
             auto& vec = is_lower ? m_lower_terms : m_upper_terms;
             if (vec.size() <= ti) {
-                vec.resize(ti + 1, constraint_bound(null_index, rational()));
+                vec.resize(ti + 1, constraint_bound(UINT_MAX, rational()));
             }
             constraint_bound& b = vec[ti];
-            if (b.first == null_index || (is_lower? b.second < v : b.second > v)) {
+            if (b.first == UINT_MAX || (is_lower? b.second < v : b.second > v)) {
                 ctx().push_trail(vector_value_trail<context, constraint_bound>(vec, ti));
                 b.first = ci;
                 b.second = v;
@@ -2045,7 +2048,7 @@ namespace smt {
                 rational val;
                 TRACE("arith", tout << vi << " " << v << "\n";);
                 if (v != null_theory_var && a.is_numeral(get_owner(v), val) && bound == val) {
-                    ci = null_constraint_index;
+                    ci = UINT_MAX;
                     return bound == val;
                 }
 
@@ -2053,7 +2056,7 @@ namespace smt {
                 if (vec.size() > ti) {
                     constraint_bound& b = vec[ti];
                     ci = b.first;
-                    return ci != null_index && bound == b.second;
+                    return ci != UINT_MAX && bound == b.second;
                 }
                 else {
                     return false;
@@ -2137,22 +2140,10 @@ namespace smt {
             if (m_solver->A_r().row_count() > m_stats.m_max_rows)
                 m_stats.m_max_rows = m_solver->A_r().row_count();
             TRACE("arith_verbose", display(tout););
-            bool print = false && m_print_counter++ % 1000 == 0;
-            stopwatch sw;
-            if (print) {
-                sw.start();
-            }
             lean::lp_status status = m_solver->find_feasible_solution();
-            if (print) {
-                sw.stop();
-            }
             m_stats.m_num_iterations = m_solver->settings().st().m_total_iterations;
             m_stats.m_num_factorizations = m_solver->settings().st().m_num_factorizations;
             m_stats.m_need_to_solve_inf = m_solver->settings().st().m_need_to_solve_inf;
-            if (print) {
-                IF_VERBOSE(0, verbose_stream() << status << " " << sw.get_seconds() << " " << m_stats.m_num_iterations << " " << m_print_counter << "\n";);
-            }
-            //m_stats.m_num_iterations_with_no_progress += m_solver->settings().st().m_iters_with_no_cost_growing;
 
             switch (status) {
             case lean::lp_status::INFEASIBLE:
@@ -2175,10 +2166,11 @@ namespace smt {
         literal_vector      m_core;
         svector<enode_pair> m_eqs;
         vector<parameter>   m_params;
-        lean::constraint_index const null_constraint_index = UINT_MAX;
+
+        // lean::constraint_index const null_constraint_index = UINT_MAX; // not sure what a correct fix is
 
         void set_evidence(lean::constraint_index idx) {
-            if (idx == null_constraint_index) {
+            if (idx == UINT_MAX) {
                 return;
             }
             switch (m_constraint_sources[idx]) {
@@ -2517,9 +2509,9 @@ namespace smt {
         }        
     };
     
-    theory_lra::theory_lra(ast_manager& m, theory_arith_params& p):
+    theory_lra::theory_lra(ast_manager& m, theory_arith_params& ap):
         theory(m.get_family_id("arith")) {
-        m_imp = alloc(imp, *this, m, p);
+        m_imp = alloc(imp, *this, m, ap);
     }    
     theory_lra::~theory_lra() {
         dealloc(m_imp);
@@ -2529,6 +2521,7 @@ namespace smt {
     }
     void theory_lra::init(context * ctx) {
         theory::init(ctx);
+        m_imp->init(ctx);
     }
     bool theory_lra::internalize_atom(app * atom, bool gate_ctx) {
         return m_imp->internalize_atom(atom, gate_ctx);

src/smt/theory_lra.h

@@ -26,8 +26,9 @@ namespace smt {
     class theory_lra : public theory, public theory_opt {
         class imp;
         imp* m_imp;
+
     public:
-        theory_lra(ast_manager& m, theory_arith_params& params);
+        theory_lra(ast_manager& m, theory_arith_params& ap);
         virtual ~theory_lra();
         virtual theory* mk_fresh(context* new_ctx);
         virtual char const* get_name() const { return "lra"; }

src/tactic/smtlogics/qflra_tactic.cpp

@@ -22,7 +22,6 @@ Notes:
 #include"solve_eqs_tactic.h"
 #include"elim_uncnstr_tactic.h"
 #include"smt_tactic.h"
-// include"mip_tactic.h"
 #include"recover_01_tactic.h"
 #include"ctx_simplify_tactic.h"
 #include"probe_arith.h"
@@ -72,5 +71,18 @@ tactic * mk_qflra_tactic(ast_manager & m, params_ref const & p) {
     //                            using_params(mk_smt_tactic(), pivot_p)),
     //                    p);
 
+#if 0
+
+    params_ref simplex_0, simplex_1, simplex_2;
+    simplex_0.set_uint("lp.simplex_strategy", 0);    
+    simplex_1.set_uint("lp.simplex_strategy", 1);    
+    simplex_2.set_uint("lp.simplex_strategy", 2);
+    
+    return par(using_params(mk_smt_tactic(), simplex_0), 
+               using_params(mk_smt_tactic(), simplex_1), 
+               using_params(mk_smt_tactic(), simplex_2));
+#else
     return using_params(using_params(mk_smt_tactic(), pivot_p), p);
+#endif
+
 }

src/test/lp.cpp

@@ -34,6 +34,10 @@ Author: Lev Nachmanson
 namespace lean {
 unsigned seed = 1;
 
+    random_gen g_rand;
+    static unsigned my_random() {
+        return g_rand();
+    }
 struct simple_column_namer:public column_namer
 {
     std::string get_column_name(unsigned j) const override {
@@ -1077,7 +1081,7 @@ bool get_double_from_args_parser(const char * option, argument_parser & args_par
 
 void update_settings(argument_parser & args_parser, lp_settings& settings) {
     unsigned n;
-    settings.m_simplex_strategy = simplex_strategy_enum::no_tableau;
+    settings.m_simplex_strategy = simplex_strategy_enum::lu;
     if (get_int_from_args_parser("--rep_frq", args_parser, n))
         settings.report_frequency = n;
     else
@@ -1104,7 +1108,7 @@ void update_settings(argument_parser & args_parser, lp_settings& settings) {
         settings.harris_feasibility_tolerance = d;
     }
     if (get_int_from_args_parser("--random_seed", args_parser, n)) {
-        settings.random_seed = n;
+        settings.random_seed(n);
     }
     if (get_int_from_args_parser("--simplex_strategy", args_parser, n)) {
         settings.simplex_strategy() = static_cast<simplex_strategy_enum>(n);

src/test/smt_reader.h

@@ -56,23 +56,27 @@ namespace lean {
         struct formula_constraint {
             lconstraint_kind m_kind;
             std::vector<std::pair<mpq, std::string>> m_coeffs;
-            mpq m_right_side = numeric_traits<mpq>::zero();
+            mpq m_right_side;
             void add_pair(mpq c, std::string name) {
                 m_coeffs.push_back(make_pair(c, name));
             }
+            formula_constraint() : m_right_side(numeric_traits<mpq>::zero()) {}
         };
 
         lisp_elem m_formula_lisp_elem;
 
         std::unordered_map<std::string, unsigned> m_name_to_var_index;
-        std::vector<formula_constraint> m_constraints;
+        std::vector<formula_constraint>           m_constraints;
-        std::string m_file_name;
+        bool                                      m_is_OK;
-        std::ifstream m_file_stream;
+        unsigned                                  m_line_number;
-        std::string m_line;
+        std::string                               m_file_name;
-        bool m_is_OK = true;
+        std::ifstream                             m_file_stream;
-        unsigned m_line_number = 0;
+        std::string                               m_line;
         smt_reader(std::string file_name):
-            m_file_name(file_name), m_file_stream(file_name) {
+            m_is_OK(true),
+            m_line_number(0),
+            m_file_name(file_name), 
+            m_file_stream(file_name) {
         }
 
         void set_error() {
@@ -364,7 +368,7 @@ namespace lean {
             if (it!= m_name_to_var_index.end())
                 return it->second;
 
-            unsigned ret= m_name_to_var_index.size();
+            unsigned ret = static_cast<unsigned>(m_name_to_var_index.size());
             m_name_to_var_index[s] = ret;
             return ret;
         }

src/util/hash.cpp

@@ -83,8 +83,8 @@ unsigned string_hash(const char * str, unsigned length, unsigned init_value) {
         Z3_fallthrough;
     case 1 : 
         a+=str[0];
-        Z3_fallthrough;
         /* case 0: nothing left to add */
+        break;
     }
     mix(a,b,c);
     /*-------------------------------------------- report the result */

src/util/lp/binary_heap_priority_queue.h

@@ -16,8 +16,7 @@ class binary_heap_priority_queue {
     // indexing for A starts from 1
     vector<unsigned> m_heap; // keeps the elements of the queue
     vector<int> m_heap_inverse; // o = m_heap[m_heap_inverse[o]]
-    unsigned m_heap_size = 0;
+    unsigned m_heap_size;
-
     // is is the child place in heap
     void swap_with_parent(unsigned i);
     void put_at(unsigned i, unsigned h);
@@ -29,7 +28,7 @@ public:
 public:
     void remove(unsigned o);
     unsigned size() const { return m_heap_size; }
-    binary_heap_priority_queue(): m_heap(1) {} // the empty constructror
+    binary_heap_priority_queue(): m_heap(1), m_heap_size(0) {} // the empty constructror
     // n is the initial queue capacity.
     // The capacity will be enlarged two times automatically if needed
     binary_heap_priority_queue(unsigned n);

src/util/lp/binary_heap_priority_queue.hpp

@@ -83,7 +83,8 @@ template <typename T> void binary_heap_priority_queue<T>::remove(unsigned o) {
 template <typename T> binary_heap_priority_queue<T>::binary_heap_priority_queue(unsigned n) :
     m_priorities(n),
     m_heap(n + 1), // because the indexing for A starts from 1
-    m_heap_inverse(n, -1)
+    m_heap_inverse(n, -1),
+    m_heap_size(0)
 { }
 
 

src/util/lp/binary_heap_upair_queue.h

@@ -20,8 +20,8 @@ template <typename  T>
 class binary_heap_upair_queue {
     binary_heap_priority_queue<T> m_q;
     std::unordered_map<upair, unsigned> m_pairs_to_index;
-    vector<upair> m_pairs; // inverse to index
+    svector<upair> m_pairs; // inverse to index
-    vector<unsigned> m_available_spots;
+    svector<unsigned> m_available_spots;
 public:
     binary_heap_upair_queue(unsigned size);
 

src/util/lp/bound_analyzer_on_row.h

@@ -18,12 +18,13 @@ namespace lean {
 class bound_analyzer_on_row {
     
     linear_combination_iterator<mpq> & m_it;
-    unsigned m_row_or_term_index;
+    bound_propagator &                 m_bp;
-    int m_column_of_u = -1; // index of an unlimited from above monoid
+    unsigned           m_row_or_term_index;
-    // -1 means that such a value is not found, -2 means that at least two of such monoids were found
+    int                m_column_of_u; // index of an unlimited from above monoid
-    int m_column_of_l = -1; // index of an unlimited from below monoid
+                               // -1 means that such a value is not found, -2 means that at least two of such monoids were found
-    impq m_rs;
+    int                m_column_of_l; // index of an unlimited from below monoid
-    bound_propagator & m_bp;
+    impq               m_rs;
+
 public :
     // constructor
     bound_analyzer_on_row(
@@ -34,9 +35,11 @@ public :
                           )
         :
         m_it(it),
+        m_bp(bp),
         m_row_or_term_index(row_or_term_index),
-        m_rs(rs),
+        m_column_of_u(-1),
-        m_bp(bp)
+        m_column_of_l(-1),
+        m_rs(rs)
     {}
 
 
@@ -250,7 +253,6 @@ public :
             if (str)
                 strict = true;
         }
-
         bound /= l_coeff;
         if (is_pos(l_coeff)) {
             limit_j(m_column_of_l, bound, true, false, strict);

src/util/lp/column_info.h

@@ -15,26 +15,26 @@ inline bool is_valid(unsigned j) { return static_cast<int>(j) >= 0;}
 template <typename T>
 class column_info {
     std::string m_name;
-    bool m_low_bound_is_set = false;
+    bool        m_low_bound_is_set;
-    bool m_low_bound_is_strict = false;
+    bool        m_low_bound_is_strict;
-    bool m_upper_bound_is_set = false;
+    bool        m_upper_bound_is_set;
-    bool m_upper_bound_is_strict = false;
+    bool        m_upper_bound_is_strict;
-    T m_low_bound;
+    T           m_low_bound;
-    T m_upper_bound;
+    T           m_upper_bound;
-    T m_cost = numeric_traits<T>::zero();
+    T           m_fixed_value;
-    T m_fixed_value;
+    bool        m_is_fixed;
-    bool m_is_fixed = false;
+    T           m_cost;
-    unsigned m_column_index = static_cast<unsigned>(-1);
+    unsigned    m_column_index;
 public:
     bool operator==(const column_info & c) const {
-        return     m_name == c.m_name &&
+        return m_name == c.m_name &&
             m_low_bound_is_set == c.m_low_bound_is_set &&
             m_low_bound_is_strict == c.m_low_bound_is_strict &&
             m_upper_bound_is_set == c.m_upper_bound_is_set&&
             m_upper_bound_is_strict == c.m_upper_bound_is_strict&&
             (!m_low_bound_is_set || m_low_bound == c.m_low_bound) &&
             (!m_upper_bound_is_set || m_upper_bound == c.m_upper_bound) &&
-            m_cost == c.m_cost&&
+            m_cost == c.m_cost &&
             m_is_fixed == c.m_is_fixed &&
             (!m_is_fixed || m_fixed_value == c.m_fixed_value) &&
             m_column_index == c.m_column_index;
@@ -44,9 +44,24 @@ public:
         m_column_index = j;
     }
     // the default constructor
-    column_info() {}
+    column_info():
+        m_low_bound_is_set(false),
+        m_low_bound_is_strict(false),
+        m_upper_bound_is_set (false),
+        m_upper_bound_is_strict (false),
+        m_is_fixed(false),
+        m_cost(numeric_traits<T>::zero()),
+        m_column_index(static_cast<unsigned>(-1))
+    {}
     
-    column_info(unsigned column_index) : m_column_index(column_index) {
+    column_info(unsigned column_index) :
+        m_low_bound_is_set(false),
+        m_low_bound_is_strict(false),
+        m_upper_bound_is_set (false),
+        m_upper_bound_is_strict (false),
+        m_is_fixed(false),
+        m_cost(numeric_traits<T>::zero()),
+        m_column_index(column_index) {
     }
 
     column_info(const column_info & ci) {

src/util/lp/conversion_helper.h

@@ -0,0 +1,43 @@
+/*
+  Copyright (c) 2013 Microsoft Corporation. All rights reserved.
+
+  Author: Lev Nachmanson
+*/
+#pragma once
+namespace lean {
+template <typename V>
+struct conversion_helper {
+    static V get_low_bound(const column_info<mpq> & ci) {
+        return V(ci.get_low_bound(), ci.low_bound_is_strict()? 1 : 0);
+    }
+
+    static V get_upper_bound(const column_info<mpq> & ci) {
+        return V(ci.get_upper_bound(), ci.upper_bound_is_strict()? -1 : 0);
+    }
+};
+
+template<>
+struct conversion_helper <double> {
+    static double get_upper_bound(const column_info<mpq> & ci) {
+        if (!ci.upper_bound_is_strict())
+            return ci.get_upper_bound().get_double();
+        double eps = 0.00001;
+        if (!ci.low_bound_is_set())
+            return ci.get_upper_bound().get_double() - eps;
+        eps = std::min((ci.get_upper_bound() - ci.get_low_bound()).get_double() / 1000, eps);
+        return ci.get_upper_bound().get_double() - eps;
+    }
+
+    static double get_low_bound(const column_info<mpq> & ci) {
+        if (!ci.low_bound_is_strict())
+            return ci.get_low_bound().get_double();
+        double eps = 0.00001;
+        if (!ci.upper_bound_is_set())
+            return ci.get_low_bound().get_double() + eps;
+        eps = std::min((ci.get_upper_bound() - ci.get_low_bound()).get_double() / 1000, eps);
+        return ci.get_low_bound().get_double() + eps;
+    }
+
+};
+
+}

src/util/lp/core_solver_pretty_printer.h

@@ -16,7 +16,6 @@ template <typename T, typename X> class lp_core_solver_base; // forward definiti
 template <typename T, typename X>
 class core_solver_pretty_printer {
     std::ostream & m_out;
-    template<typename A> using vector = vector<A>;
     typedef std::string string;
     lp_core_solver_base<T, X> & m_core_solver;
     vector<unsigned> m_column_widths;
@@ -34,15 +33,15 @@ class core_solver_pretty_printer {
     std::string m_cost_title;
     std::string m_basis_heading_title;
     std::string m_x_title;
-    std::string m_low_bounds_title = "low";
+    std::string m_low_bounds_title;
-    std::string m_upp_bounds_title = "upp";
+    std::string m_upp_bounds_title;
-    std::string m_exact_norm_title = "exact cn";
+    std::string m_exact_norm_title;
-    std::string m_approx_norm_title = "approx cn";
+    std::string m_approx_norm_title;
 
 
     unsigned ncols() { return m_core_solver.m_A.column_count(); }
     unsigned nrows() { return m_core_solver.m_A.row_count(); }
-    unsigned m_artificial_start = std::numeric_limits<unsigned>::max();
+    unsigned m_artificial_start;
     indexed_vector<T> m_w_buff;
     indexed_vector<T> m_ed_buff;
     vector<T> m_exact_column_norms;

src/util/lp/core_solver_pretty_printer.hpp

@@ -23,6 +23,12 @@ core_solver_pretty_printer<T, X>::core_solver_pretty_printer(lp_core_solver_base
     m_rs(ncols(), zero_of_type<X>()),
     m_w_buff(core_solver.m_w),
     m_ed_buff(core_solver.m_ed) {
+    m_low_bounds_title = "low";
+    m_upp_bounds_title = "upp";
+    m_exact_norm_title = "exact cn";
+    m_approx_norm_title = "approx cn";
+    m_artificial_start = std::numeric_limits<unsigned>::max();
+
     m_column_widths.resize(core_solver.m_A.column_count(), 0),
     init_m_A_and_signs();
     init_costs();

src/util/lp/init_lar_solver.h

@@ -0,0 +1,576 @@
+/*
+  Copyright (c) 2017 Microsoft Corporation
+  Author: Lev Nachmanson
+*/
+
+// here we are inside lean::lar_solver class
+
+bool strategy_is_undecided() const {
+    return m_settings.simplex_strategy() == simplex_strategy_enum::undecided;
+}
+
+var_index add_var(unsigned ext_j) {
+    var_index i;
+    lean_assert (ext_j < m_terms_start_index); 
+
+    if (ext_j >= m_terms_start_index)
+        throw 0; // todo : what is the right way to exit?
+            
+    if (try_get_val(m_ext_vars_to_columns, ext_j, i)) {
+        return i;
+    }
+    lean_assert(m_vars_to_ul_pairs.size() == A_r().column_count());
+    i = A_r().column_count();
+    m_vars_to_ul_pairs.push_back (ul_pair(static_cast<unsigned>(-1)));
+    add_non_basic_var_to_core_fields(ext_j);
+    lean_assert(sizes_are_correct());
+    return i;
+}
+
+void register_new_ext_var_index(unsigned ext_v) {
+    lean_assert(!contains(m_ext_vars_to_columns, ext_v));
+    unsigned j = static_cast<unsigned>(m_ext_vars_to_columns.size());
+    m_ext_vars_to_columns[ext_v] = j;
+    lean_assert(m_columns_to_ext_vars_or_term_indices.size() == j);
+    m_columns_to_ext_vars_or_term_indices.push_back(ext_v);
+}
+
+void add_non_basic_var_to_core_fields(unsigned ext_j) {
+    register_new_ext_var_index(ext_j);
+    m_mpq_lar_core_solver.m_column_types.push_back(column_type::free_column);
+    m_columns_with_changed_bound.increase_size_by_one();
+    add_new_var_to_core_fields_for_mpq(false);
+    if (use_lu())
+        add_new_var_to_core_fields_for_doubles(false);
+}
+
+void add_new_var_to_core_fields_for_doubles(bool register_in_basis) {
+    unsigned j = A_d().column_count();
+    A_d().add_column();
+    lean_assert(m_mpq_lar_core_solver.m_d_x.size() == j);
+    //        lean_assert(m_mpq_lar_core_solver.m_d_low_bounds.size() == j && m_mpq_lar_core_solver.m_d_upper_bounds.size() == j);  // restore later
+    m_mpq_lar_core_solver.m_d_x.resize(j + 1 ); 
+    m_mpq_lar_core_solver.m_d_low_bounds.resize(j + 1);
+    m_mpq_lar_core_solver.m_d_upper_bounds.resize(j + 1);
+    lean_assert(m_mpq_lar_core_solver.m_d_heading.size() == j); // as A().column_count() on the entry to the method
+    if (register_in_basis) {
+        A_d().add_row();
+        m_mpq_lar_core_solver.m_d_heading.push_back(m_mpq_lar_core_solver.m_d_basis.size());
+        m_mpq_lar_core_solver.m_d_basis.push_back(j);
+    }else {
+        m_mpq_lar_core_solver.m_d_heading.push_back(- static_cast<int>(m_mpq_lar_core_solver.m_d_nbasis.size()) - 1);
+        m_mpq_lar_core_solver.m_d_nbasis.push_back(j);
+    }
+}
+
+void add_new_var_to_core_fields_for_mpq(bool register_in_basis) {
+    unsigned j = A_r().column_count();
+    A_r().add_column();
+    lean_assert(m_mpq_lar_core_solver.m_r_x.size() == j);
+    //        lean_assert(m_mpq_lar_core_solver.m_r_low_bounds.size() == j && m_mpq_lar_core_solver.m_r_upper_bounds.size() == j);  // restore later
+    m_mpq_lar_core_solver.m_r_x.resize(j + 1);
+    m_mpq_lar_core_solver.m_r_low_bounds.increase_size_by_one();
+    m_mpq_lar_core_solver.m_r_upper_bounds.increase_size_by_one();
+    m_mpq_lar_core_solver.m_r_solver.m_inf_set.increase_size_by_one();
+    m_mpq_lar_core_solver.m_r_solver.m_costs.resize(j + 1);
+    m_mpq_lar_core_solver.m_r_solver.m_d.resize(j + 1);
+    lean_assert(m_mpq_lar_core_solver.m_r_heading.size() == j); // as A().column_count() on the entry to the method
+    if (register_in_basis) {
+        A_r().add_row();
+        m_mpq_lar_core_solver.m_r_heading.push_back(m_mpq_lar_core_solver.m_r_basis.size());
+        m_mpq_lar_core_solver.m_r_basis.push_back(j);
+        if (m_settings.bound_propagation())
+            m_rows_with_changed_bounds.insert(A_r().row_count() - 1);
+    } else {
+        m_mpq_lar_core_solver.m_r_heading.push_back(- static_cast<int>(m_mpq_lar_core_solver.m_r_nbasis.size()) - 1);
+        m_mpq_lar_core_solver.m_r_nbasis.push_back(j);
+    }
+}
+
+
+var_index add_term_undecided(const vector<std::pair<mpq, var_index>> & coeffs,
+                             const mpq &m_v) {
+    m_terms.push_back(new lar_term(coeffs, m_v));
+    m_orig_terms.push_back(new lar_term(coeffs, m_v));
+    return m_terms_start_index +  m_terms.size() - 1;
+}
+
+// terms
+var_index add_term(const vector<std::pair<mpq, var_index>> & coeffs,
+                   const mpq &m_v) {
+    if (strategy_is_undecided())
+        return add_term_undecided(coeffs, m_v);
+
+    m_terms.push_back(new lar_term(coeffs, m_v));
+    m_orig_terms.push_back(new lar_term(coeffs, m_v));
+    unsigned adjusted_term_index = m_terms.size() - 1;
+    var_index ret = m_terms_start_index + adjusted_term_index;
+    if (use_tableau() && !coeffs.empty()) {
+        add_row_for_term(m_orig_terms.back(), ret);
+        if (m_settings.bound_propagation())
+            m_rows_with_changed_bounds.insert(A_r().row_count() - 1);
+    } 
+    lean_assert(m_ext_vars_to_columns.size() == A_r().column_count());
+    return ret;
+}
+
+void add_row_for_term(const lar_term * term, unsigned term_ext_index) {
+    lean_assert(sizes_are_correct());
+    add_row_from_term_no_constraint(term, term_ext_index);
+    lean_assert(sizes_are_correct());
+}
+
+void add_row_from_term_no_constraint(const lar_term * term, unsigned term_ext_index) {
+    register_new_ext_var_index(term_ext_index);
+    // j will be a new variable
+	unsigned j = A_r().column_count();
+    ul_pair ul(j);
+    m_vars_to_ul_pairs.push_back(ul);
+    add_basic_var_to_core_fields();
+    if (use_tableau()) {
+        auto it = iterator_on_term_with_basis_var(*term, j);
+        A_r().fill_last_row_with_pivoting(it,
+                                          m_mpq_lar_core_solver.m_r_solver.m_basis_heading);
+        m_mpq_lar_core_solver.m_r_solver.m_b.resize(A_r().column_count(), zero_of_type<mpq>());
+    } else {
+        fill_last_row_of_A_r(A_r(), term);
+    }
+    m_mpq_lar_core_solver.m_r_x[j] = get_basic_var_value_from_row_directly(A_r().row_count() - 1);
+    if (use_lu())
+        fill_last_row_of_A_d(A_d(), term);
+}
+
+void add_basic_var_to_core_fields() {
+    bool use_lu = m_mpq_lar_core_solver.need_to_presolve_with_double_solver();
+    lean_assert(!use_lu || A_r().column_count() == A_d().column_count());
+    m_mpq_lar_core_solver.m_column_types.push_back(column_type::free_column);
+    m_columns_with_changed_bound.increase_size_by_one();
+    m_rows_with_changed_bounds.increase_size_by_one();
+    add_new_var_to_core_fields_for_mpq(true);
+    if (use_lu)
+        add_new_var_to_core_fields_for_doubles(true);
+}
+
+constraint_index add_var_bound(var_index j, lconstraint_kind kind, const mpq & right_side)  {
+	constraint_index ci = m_constraints.size();
+    if (!is_term(j)) { // j is a var
+        auto vc = new lar_var_constraint(j, kind, right_side);
+        m_constraints.push_back(vc);
+        update_column_type_and_bound(j, kind, right_side, ci);
+    } else {
+        add_var_bound_on_constraint_for_term(j, kind, right_side, ci);
+    }
+    lean_assert(sizes_are_correct());
+    return ci;
+}
+
+void update_column_type_and_bound(var_index j, lconstraint_kind kind, const mpq & right_side, constraint_index constr_index) {
+    switch(m_mpq_lar_core_solver.m_column_types[j]) {
+    case column_type::free_column:
+        update_free_column_type_and_bound(j, kind, right_side, constr_index);
+        break;
+    case column_type::boxed:
+        update_boxed_column_type_and_bound(j, kind, right_side, constr_index);
+        break;
+    case column_type::low_bound:
+        update_low_bound_column_type_and_bound(j, kind, right_side, constr_index);
+        break;
+    case column_type::upper_bound:
+        update_upper_bound_column_type_and_bound(j, kind, right_side, constr_index);
+        break;
+    case column_type::fixed:
+        update_fixed_column_type_and_bound(j, kind, right_side, constr_index);
+        break;
+    default:
+        lean_assert(false); // cannot be here
+    }
+}
+
+void add_var_bound_on_constraint_for_term(var_index j, lconstraint_kind kind, const mpq & right_side, constraint_index ci) {
+    lean_assert(is_term(j));
+    unsigned adjusted_term_index = adjust_term_index(j);
+    unsigned term_j;
+    if (try_get_val(m_ext_vars_to_columns, j, term_j)) {
+        mpq rs = right_side - m_orig_terms[adjusted_term_index]->m_v;
+        m_constraints.push_back(new lar_term_constraint(m_orig_terms[adjusted_term_index], kind, right_side));
+        update_column_type_and_bound(term_j, kind, rs, ci);
+    }
+    else {
+        add_constraint_from_term_and_create_new_column_row(j, m_orig_terms[adjusted_term_index], kind, right_side);
+    }
+}
+
+
+void add_constraint_from_term_and_create_new_column_row(unsigned term_j, const lar_term* term,
+                                                        lconstraint_kind kind, const mpq & right_side) {
+
+    add_row_from_term_no_constraint(term, term_j);
+    unsigned j = A_r().column_count() - 1;
+    update_column_type_and_bound(j, kind, right_side - term->m_v, m_constraints.size());
+    m_constraints.push_back(new lar_term_constraint(term, kind, right_side));
+    lean_assert(A_r().column_count() == m_mpq_lar_core_solver.m_r_solver.m_costs.size());
+}
+
+void decide_on_strategy_and_adjust_initial_state() {
+	lean_assert(strategy_is_undecided());
+	if (m_vars_to_ul_pairs.size() > m_settings.column_number_threshold_for_using_lu_in_lar_solver) {
+        m_settings.simplex_strategy() = simplex_strategy_enum::lu;
+    } else {
+        m_settings.simplex_strategy() = simplex_strategy_enum::tableau_rows; // todo: when to switch to tableau_costs?
+    }
+    adjust_initial_state();
+}
+
+void adjust_initial_state() {
+    switch (m_settings.simplex_strategy()) {
+    case simplex_strategy_enum::lu:
+        adjust_initial_state_for_lu();
+        break;
+    case simplex_strategy_enum::tableau_rows:
+        adjust_initial_state_for_tableau_rows();
+        break;
+    case simplex_strategy_enum::tableau_costs:
+        lean_assert(false); // not implemented
+    case simplex_strategy_enum::undecided:
+        adjust_initial_state_for_tableau_rows();
+        break;
+    }
+}
+
+void adjust_initial_state_for_lu() {
+    copy_from_mpq_matrix(A_d());
+	unsigned n = A_d().column_count();
+	m_mpq_lar_core_solver.m_d_x.resize(n);
+	m_mpq_lar_core_solver.m_d_low_bounds.resize(n);
+	m_mpq_lar_core_solver.m_d_upper_bounds.resize(n);
+	m_mpq_lar_core_solver.m_d_heading = m_mpq_lar_core_solver.m_r_heading;
+	m_mpq_lar_core_solver.m_d_basis = m_mpq_lar_core_solver.m_r_basis;
+
+	/*
+    unsigned j = A_d().column_count();
+    A_d().add_column();
+    lean_assert(m_mpq_lar_core_solver.m_d_x.size() == j);
+    //        lean_assert(m_mpq_lar_core_solver.m_d_low_bounds.size() == j && m_mpq_lar_core_solver.m_d_upper_bounds.size() == j);  // restore later
+    m_mpq_lar_core_solver.m_d_x.resize(j + 1 ); 
+    m_mpq_lar_core_solver.m_d_low_bounds.resize(j + 1);
+    m_mpq_lar_core_solver.m_d_upper_bounds.resize(j + 1);
+    lean_assert(m_mpq_lar_core_solver.m_d_heading.size() == j); // as A().column_count() on the entry to the method
+    if (register_in_basis) {
+        A_d().add_row();
+        m_mpq_lar_core_solver.m_d_heading.push_back(m_mpq_lar_core_solver.m_d_basis.size());
+        m_mpq_lar_core_solver.m_d_basis.push_back(j);
+    }else {
+        m_mpq_lar_core_solver.m_d_heading.push_back(- static_cast<int>(m_mpq_lar_core_solver.m_d_nbasis.size()) - 1);
+        m_mpq_lar_core_solver.m_d_nbasis.push_back(j);
+        }*/
+}
+
+void adjust_initial_state_for_tableau_rows() {
+    for (unsigned j = 0; j < m_terms.size(); j++) {
+        if (contains(m_ext_vars_to_columns, j + m_terms_start_index))
+            continue;
+        add_row_from_term_no_constraint(m_terms[j], j + m_terms_start_index);
+    }
+}
+
+// this fills the last row of A_d and sets the basis column: -1 in the last column of the row
+void fill_last_row_of_A_d(static_matrix<double, double> & A, const lar_term* ls) {
+    lean_assert(A.row_count() > 0);
+    lean_assert(A.column_count() > 0);
+    unsigned last_row = A.row_count() - 1;
+    lean_assert(A.m_rows[last_row].empty());
+
+    for (auto & t : ls->m_coeffs) {
+        lean_assert(!is_zero(t.second));
+        var_index j = t.first;
+        A.set(last_row, j, - t.second.get_double());
+    }
+
+    unsigned basis_j = A.column_count() - 1;
+    A.set(last_row, basis_j, - 1 );
+}
+
+void update_free_column_type_and_bound(var_index j, lconstraint_kind kind, const mpq & right_side, constraint_index constr_ind) {
+    mpq y_of_bound(0);
+    switch (kind) {
+    case LT:
+        y_of_bound = -1;
+    case LE:
+        m_mpq_lar_core_solver.m_column_types[j] = column_type::upper_bound;
+        lean_assert(m_mpq_lar_core_solver.m_column_types()[j] == column_type::upper_bound);
+        lean_assert(m_mpq_lar_core_solver.m_r_upper_bounds.size() > j);
+        {
+            auto up = numeric_pair<mpq>(right_side, y_of_bound);
+            m_mpq_lar_core_solver.m_r_upper_bounds[j] = up;
+        }
+        set_upper_bound_witness(j, constr_ind);
+        break;
+    case GT:
+        y_of_bound = 1;
+    case GE:
+        m_mpq_lar_core_solver.m_column_types[j] = column_type::low_bound;
+        lean_assert(m_mpq_lar_core_solver.m_r_upper_bounds.size() > j);
+        {
+            auto low = numeric_pair<mpq>(right_side, y_of_bound);
+            m_mpq_lar_core_solver.m_r_low_bounds[j] = low;
+        }
+        set_low_bound_witness(j, constr_ind);
+        break;
+    case EQ:
+        m_mpq_lar_core_solver.m_column_types[j] = column_type::fixed;
+        m_mpq_lar_core_solver.m_r_low_bounds[j] = m_mpq_lar_core_solver.m_r_upper_bounds[j] = numeric_pair<mpq>(right_side, zero_of_type<mpq>());
+        set_upper_bound_witness(j, constr_ind);
+        set_low_bound_witness(j, constr_ind);
+        break;
+
+    default:
+        lean_unreachable();
+                
+    }
+    m_columns_with_changed_bound.insert(j);
+}
+
+void update_upper_bound_column_type_and_bound(var_index j, lconstraint_kind kind, const mpq & right_side, constraint_index ci) {
+    lean_assert(m_mpq_lar_core_solver.m_column_types()[j] == column_type::upper_bound);
+    mpq y_of_bound(0);
+    switch (kind) {
+    case LT:
+        y_of_bound = -1;
+    case LE:
+        {
+            auto up = numeric_pair<mpq>(right_side, y_of_bound);
+            if (up < m_mpq_lar_core_solver.m_r_upper_bounds()[j]) {
+                m_mpq_lar_core_solver.m_r_upper_bounds[j] = up;
+                set_upper_bound_witness(j, ci);
+                m_columns_with_changed_bound.insert(j);
+            }
+        }
+        break;
+    case GT:
+        y_of_bound = 1;
+    case GE:            
+        m_mpq_lar_core_solver.m_column_types[j] = column_type::boxed;
+        {
+            auto low = numeric_pair<mpq>(right_side, y_of_bound);
+            m_mpq_lar_core_solver.m_r_low_bounds[j] = low;
+            set_low_bound_witness(j, ci);
+            m_columns_with_changed_bound.insert(j);
+            if (low > m_mpq_lar_core_solver.m_r_upper_bounds[j]) {
+                m_status = INFEASIBLE;
+                m_infeasible_column_index = j;
+            } else {
+                m_mpq_lar_core_solver.m_column_types[j] = m_mpq_lar_core_solver.m_r_low_bounds()[j] < m_mpq_lar_core_solver.m_r_upper_bounds()[j]? column_type::boxed : column_type::fixed;
+            }                     
+        }
+        break;
+    case EQ:
+        {
+            auto v = numeric_pair<mpq>(right_side, zero_of_type<mpq>());
+            if (v > m_mpq_lar_core_solver.m_r_upper_bounds[j]) {
+                m_status = INFEASIBLE;
+                set_low_bound_witness(j, ci);
+                m_infeasible_column_index = j;
+            } else {
+                m_mpq_lar_core_solver.m_r_low_bounds[j] = m_mpq_lar_core_solver.m_r_upper_bounds[j] = v;
+                m_columns_with_changed_bound.insert(j);
+                set_low_bound_witness(j, ci);
+                set_upper_bound_witness(j, ci);
+                m_mpq_lar_core_solver.m_column_types[j] = column_type::fixed;
+            }
+            break;
+        }
+        break;
+
+    default:
+        lean_unreachable();
+                
+    }
+}
+    
+void update_boxed_column_type_and_bound(var_index j, lconstraint_kind kind, const mpq & right_side, constraint_index ci) {
+    lean_assert(m_status == INFEASIBLE || (m_mpq_lar_core_solver.m_column_types()[j] == column_type::boxed && m_mpq_lar_core_solver.m_r_low_bounds()[j] < m_mpq_lar_core_solver.m_r_upper_bounds()[j]));
+    mpq y_of_bound(0);
+    switch (kind) {
+    case LT:
+        y_of_bound = -1;
+    case LE:
+        {
+            auto up = numeric_pair<mpq>(right_side, y_of_bound);
+            if (up < m_mpq_lar_core_solver.m_r_upper_bounds[j]) {
+                m_mpq_lar_core_solver.m_r_upper_bounds[j] = up;
+                set_upper_bound_witness(j, ci);
+                m_columns_with_changed_bound.insert(j);
+            }
+
+            if (up < m_mpq_lar_core_solver.m_r_low_bounds[j]) {
+                m_status = INFEASIBLE;
+                lean_assert(false);
+                m_infeasible_column_index = j;
+            } else {
+                if (m_mpq_lar_core_solver.m_r_low_bounds()[j] == m_mpq_lar_core_solver.m_r_upper_bounds()[j])
+                    m_mpq_lar_core_solver.m_column_types[j] = column_type::fixed;
+            }                    
+        }
+        break;
+    case GT:
+        y_of_bound = 1;
+    case GE:            
+        {
+            auto low = numeric_pair<mpq>(right_side, y_of_bound);
+            if (low > m_mpq_lar_core_solver.m_r_low_bounds[j]) {
+                m_mpq_lar_core_solver.m_r_low_bounds[j] = low;
+                m_columns_with_changed_bound.insert(j);
+                set_low_bound_witness(j, ci);
+            }
+            if (low > m_mpq_lar_core_solver.m_r_upper_bounds[j]) {
+                m_status = INFEASIBLE;
+                m_infeasible_column_index = j;
+            } else if ( low == m_mpq_lar_core_solver.m_r_upper_bounds[j]) {
+                m_mpq_lar_core_solver.m_column_types[j] = column_type::fixed;
+            }
+        }
+        break;
+    case EQ:
+        {
+            auto v = numeric_pair<mpq>(right_side, zero_of_type<mpq>());
+            if (v < m_mpq_lar_core_solver.m_r_low_bounds[j]) {
+                m_status = INFEASIBLE;
+                m_infeasible_column_index = j;
+                set_upper_bound_witness(j, ci);                    
+            } else if (v > m_mpq_lar_core_solver.m_r_upper_bounds[j]) {
+                m_status = INFEASIBLE;
+                m_infeasible_column_index = j;
+                set_low_bound_witness(j, ci);                    
+            } else {
+                m_mpq_lar_core_solver.m_r_low_bounds[j] = m_mpq_lar_core_solver.m_r_upper_bounds[j] = v;
+                set_low_bound_witness(j, ci);
+                set_upper_bound_witness(j, ci);
+                m_mpq_lar_core_solver.m_column_types[j] = column_type::fixed;
+                m_columns_with_changed_bound.insert(j);
+            }
+                
+            break;
+        }
+
+    default:
+        lean_unreachable();
+                
+    }
+}
+void update_low_bound_column_type_and_bound(var_index j, lconstraint_kind kind, const mpq & right_side, constraint_index ci) {
+    lean_assert(m_mpq_lar_core_solver.m_column_types()[j] == column_type::low_bound);
+    mpq y_of_bound(0);
+    switch (kind) {
+    case LT:
+        y_of_bound = -1;
+    case LE:
+        {
+            auto up = numeric_pair<mpq>(right_side, y_of_bound);
+            m_mpq_lar_core_solver.m_r_upper_bounds[j] = up;
+            set_upper_bound_witness(j, ci);
+            m_columns_with_changed_bound.insert(j);
+
+            if (up < m_mpq_lar_core_solver.m_r_low_bounds[j]) {
+                m_status = INFEASIBLE;
+                m_infeasible_column_index = j;
+            } else {
+                m_mpq_lar_core_solver.m_column_types[j] = m_mpq_lar_core_solver.m_r_low_bounds()[j] < m_mpq_lar_core_solver.m_r_upper_bounds()[j]? column_type::boxed : column_type::fixed;
+            }                    
+        }
+        break;
+    case GT:
+        y_of_bound = 1;
+    case GE:            
+        {
+            auto low = numeric_pair<mpq>(right_side, y_of_bound);
+            if (low > m_mpq_lar_core_solver.m_r_low_bounds[j]) {
+                m_mpq_lar_core_solver.m_r_low_bounds[j] = low;
+                m_columns_with_changed_bound.insert(j);
+                set_low_bound_witness(j, ci);
+            }
+        }
+        break;
+    case EQ:
+        {
+            auto v = numeric_pair<mpq>(right_side, zero_of_type<mpq>());
+            if (v < m_mpq_lar_core_solver.m_r_low_bounds[j]) {
+                m_status = INFEASIBLE;
+                m_infeasible_column_index = j;
+                set_upper_bound_witness(j, ci);                    
+            } else {
+                m_mpq_lar_core_solver.m_r_low_bounds[j] = m_mpq_lar_core_solver.m_r_upper_bounds[j] = v;
+                set_low_bound_witness(j, ci);
+                set_upper_bound_witness(j, ci);
+                m_mpq_lar_core_solver.m_column_types[j] = column_type::fixed;
+            }
+            m_columns_with_changed_bound.insert(j);
+            break;
+        }
+
+    default:
+        lean_unreachable();
+                
+    }
+}
+
+void update_fixed_column_type_and_bound(var_index j, lconstraint_kind kind, const mpq & right_side, constraint_index ci) {
+    lean_assert(m_status == INFEASIBLE || (m_mpq_lar_core_solver.m_column_types()[j] == column_type::fixed && m_mpq_lar_core_solver.m_r_low_bounds()[j] == m_mpq_lar_core_solver.m_r_upper_bounds()[j]));
+    lean_assert(m_status == INFEASIBLE || (m_mpq_lar_core_solver.m_r_low_bounds()[j].y.is_zero() && m_mpq_lar_core_solver.m_r_upper_bounds()[j].y.is_zero()));
+    auto v = numeric_pair<mpq>(right_side, mpq(0));
+        
+    mpq y_of_bound(0);
+    switch (kind) {
+    case LT:
+        if (v <= m_mpq_lar_core_solver.m_r_low_bounds[j]) {
+            m_status = INFEASIBLE;
+            m_infeasible_column_index = j;
+            set_upper_bound_witness(j, ci);
+        }                   
+        break;
+    case LE:
+        {
+            if (v < m_mpq_lar_core_solver.m_r_low_bounds[j]) {
+                m_status = INFEASIBLE;
+                m_infeasible_column_index = j;
+                set_upper_bound_witness(j, ci);
+            }                   
+        }
+        break;
+    case GT:
+        {
+            if (v >= m_mpq_lar_core_solver.m_r_upper_bounds[j]) {
+                m_status = INFEASIBLE;
+                m_infeasible_column_index =j;
+                set_low_bound_witness(j, ci);
+            }
+        }
+        break;
+    case GE:            
+        {
+            if (v > m_mpq_lar_core_solver.m_r_upper_bounds[j]) {
+                m_status = INFEASIBLE;
+                m_infeasible_column_index = j;
+                set_low_bound_witness(j, ci);
+            }
+        }
+        break;
+    case EQ:
+        {
+            if (v < m_mpq_lar_core_solver.m_r_low_bounds[j]) {
+                m_status = INFEASIBLE;
+                m_infeasible_column_index = j;
+                set_upper_bound_witness(j, ci);                    
+            } else if (v > m_mpq_lar_core_solver.m_r_upper_bounds[j]) {
+                m_status = INFEASIBLE;
+                m_infeasible_column_index = j;
+                set_low_bound_witness(j, ci);                    
+            } 
+            break;
+        }
+
+    default:
+        lean_unreachable();
+                
+    }
+}
+    

src/util/lp/iterator_on_column.h

@@ -11,7 +11,7 @@ template <typename T, typename X>
 struct iterator_on_column:linear_combination_iterator<T> {
     const vector<column_cell>& m_column; // the offset in term coeffs
     const static_matrix<T, X> & m_A;
-    int m_i = -1; // the initial offset in the column
+    int m_i; // the initial offset in the column
     unsigned size() const { return m_column.size(); }
     iterator_on_column(const vector<column_cell>& column, const static_matrix<T,X> & A) // the offset in term coeffs
         :

src/util/lp/iterator_on_indexed_vector.h

@@ -8,8 +8,11 @@ namespace lean {
 template <typename T>
 struct iterator_on_indexed_vector:linear_combination_iterator<T> {
     const indexed_vector<T> & m_v;
-    unsigned m_offset = 0;
+    unsigned m_offset;
-    iterator_on_indexed_vector(const indexed_vector<T> & v) : m_v(v){}
+    iterator_on_indexed_vector(const indexed_vector<T> & v) :
+        m_v(v),
+        m_offset(0)
+    {}
     unsigned size() const { return m_v.m_index.size(); }
     bool next(T & a, unsigned & i) {
         if (m_offset >= m_v.m_index.size())

src/util/lp/iterator_on_pivot_row.h

@@ -7,12 +7,14 @@
 namespace lean {
 template <typename T>
 struct iterator_on_pivot_row:linear_combination_iterator<T> {
-    bool m_basis_returned = false;
+    bool m_basis_returned;
     const indexed_vector<T> & m_v;
     unsigned m_basis_j;
     iterator_on_indexed_vector<T> m_it;
     unsigned size() const { return m_it.size(); }
-    iterator_on_pivot_row(const indexed_vector<T> & v, unsigned basis_j) : m_v(v), m_basis_j(basis_j), m_it(v) {}
+    iterator_on_pivot_row(const indexed_vector<T> & v, unsigned basis_j) :
+        m_basis_returned(false),
+        m_v(v), m_basis_j(basis_j), m_it(v) {}
     bool next(T & a, unsigned & i) {
         if (m_basis_returned == false) {
             m_basis_returned = true;

src/util/lp/iterator_on_row.h

@@ -8,8 +8,8 @@ namespace lean {
 template <typename T>
 struct iterator_on_row:linear_combination_iterator<T> {
     const vector<row_cell<T>> & m_row;
-    unsigned m_i= 0; // offset
+    unsigned m_i; // offset
-    iterator_on_row(const vector<row_cell<T>> & row) : m_row(row)
+    iterator_on_row(const vector<row_cell<T>> & row) : m_row(row), m_i(0)
     {}
     unsigned size() const { return m_row.size(); }
     bool next(T & a, unsigned & i) {

src/util/lp/iterator_on_term_with_basis_var.h

@@ -8,14 +8,15 @@
 #include "util/lp/lar_term.h"
 namespace lean {
 struct iterator_on_term_with_basis_var:linear_combination_iterator<mpq> {
-    std::unordered_map<unsigned, mpq>::const_iterator m_i; // the offset in term coeffs
-    bool m_term_j_returned = false;
     const lar_term & m_term;
-    unsigned m_term_j;
+    std::unordered_map<unsigned, mpq>::const_iterator m_i; // the offset in term coeffs
+    bool             m_term_j_returned;
+    unsigned         m_term_j;
     unsigned size() const {return static_cast<unsigned>(m_term.m_coeffs.size() + 1);}
     iterator_on_term_with_basis_var(const lar_term & t, unsigned term_j) :
-        m_i(t.m_coeffs.begin()),
         m_term(t),
+        m_i(t.m_coeffs.begin()),
+        m_term_j_returned(false),
         m_term_j(term_j) {}
 
     bool next(mpq & a, unsigned & i) {

src/util/lp/lar_core_solver.h

@@ -25,7 +25,7 @@ class lar_core_solver  {
     // to grow and is set to -1  otherwise
     int m_sign_of_entering_delta;
     vector<std::pair<mpq, unsigned>> m_infeasible_linear_combination;
-    int m_infeasible_sum_sign = 0; // todo: get rid of this field
+    int m_infeasible_sum_sign; // todo: get rid of this field
     vector<numeric_pair<mpq>> m_right_sides_dummy;
     vector<mpq> m_costs_dummy;
     vector<double> m_d_right_sides_dummy;
@@ -216,8 +216,6 @@ public:
 
     
     void pop(unsigned k) {
-        m_stacked_simplex_strategy.pop(k);
-        bool use_tableau = m_stacked_simplex_strategy() != simplex_strategy_enum::no_tableau;
         // rationals
         if (!settings().use_tableau()) 
             m_r_A.pop(k);
@@ -232,7 +230,7 @@ public:
         m_r_x.resize(m_r_A.column_count());
         m_r_solver.m_costs.resize(m_r_A.column_count());
         m_r_solver.m_d.resize(m_r_A.column_count());
-        if(!use_tableau)
+        if(!settings().use_tableau())
             pop_markowitz_counts(k);
         m_d_A.pop(k);
         if (m_d_solver.m_factorization != nullptr) {
@@ -242,13 +240,14 @@ public:
         
         m_d_x.resize(m_d_A.column_count());
         pop_basis(k);
-
+        m_stacked_simplex_strategy.pop(k);
+        settings().simplex_strategy() = m_stacked_simplex_strategy;
         lean_assert(m_r_solver.basis_heading_is_correct());
         lean_assert(!need_to_presolve_with_double_solver() || m_d_solver.basis_heading_is_correct());
     }
 
     bool need_to_presolve_with_double_solver() const {
-        return settings().presolve_with_double_solver_for_lar && !settings().use_tableau();
+        return settings().simplex_strategy() == simplex_strategy_enum::lu;
     }
 
     template <typename L>
@@ -368,7 +367,7 @@ public:
                       s.m_x[j] = s.m_low_bounds[j];
                       break;
                   case column_type::boxed:
-                      if (my_random() % 2) {
+                      if (settings().random_next() % 2) {
                           s.m_x[j] = s.m_low_bounds[j];
                       } else {
                           s.m_x[j] = s.m_upper_bounds[j];
@@ -600,7 +599,7 @@ public:
         }
 
         if (no_r_lu()) { // it is the case where m_d_solver gives a degenerated basis, we need to roll back
-            std::cout << "no_r_lu" << std::endl;
+            //            std::cout << "no_r_lu" << std::endl;
             catch_up_in_lu_in_reverse(changes_of_basis, m_r_solver);
             m_r_solver.find_feasible_solution();
             m_d_basis = m_r_basis;
@@ -774,8 +773,8 @@ public:
     }
 
     
-    mpq find_delta_for_strict_bounds() const{
+    mpq find_delta_for_strict_bounds(const mpq & initial_delta) const{
-        mpq delta = numeric_traits<mpq>::one();
+        mpq delta = initial_delta;
         for (unsigned j = 0; j < m_r_A.column_count(); j++ ) {
             if (low_bound_is_set(j))
                 update_delta(delta, m_r_low_bounds[j], m_r_x[j]);

src/util/lp/lar_core_solver.hpp

@@ -14,7 +14,8 @@ namespace lean {
 lar_core_solver::lar_core_solver(
                                  lp_settings & settings,
                                  const column_namer & column_names
-):
+                                 ):
+    m_infeasible_sum_sign(0),
     m_r_solver(m_r_A,
                     m_right_sides_dummy,
                     m_r_x,

src/util/lp/lar_solver.h

@@ -29,86 +29,38 @@
 #include "util/lp/iterator_on_term_with_basis_var.h"
 #include "util/lp/iterator_on_row.h"
 #include "util/lp/quick_xplain.h"
+#include "util/lp/conversion_helper.h"
 namespace lean {
-template <typename V>
-struct conversion_helper {
-    static V get_low_bound(const column_info<mpq> & ci) {
-        return V(ci.get_low_bound(), ci.low_bound_is_strict()? 1 : 0);
-    }
-
-    static V get_upper_bound(const column_info<mpq> & ci) {
-        return V(ci.get_upper_bound(), ci.upper_bound_is_strict()? -1 : 0);
-    }
-};
-
-template<>
-struct conversion_helper <double> {
-    static double get_upper_bound(const column_info<mpq> & ci) {
-        if (!ci.upper_bound_is_strict())
-            return ci.get_upper_bound().get_double();
-        double eps = 0.00001;
-        if (!ci.low_bound_is_set())
-            return ci.get_upper_bound().get_double() - eps;
-        eps = std::min((ci.get_upper_bound() - ci.get_low_bound()).get_double() / 1000, eps);
-        return ci.get_upper_bound().get_double() - eps;
-    }
-
-    static double get_low_bound(const column_info<mpq> & ci) {
-        if (!ci.low_bound_is_strict())
-            return ci.get_low_bound().get_double();
-        double eps = 0.00001;
-        if (!ci.upper_bound_is_set())
-            return ci.get_low_bound().get_double() + eps;
-        eps = std::min((ci.get_upper_bound() - ci.get_low_bound()).get_double() / 1000, eps);
-        return ci.get_low_bound().get_double() + eps;
-    }
-
-};
-
-struct constraint_index_and_column_struct {
-    int m_ci = -1;
-    int m_j = -1;
-    constraint_index_and_column_struct() {}
-    constraint_index_and_column_struct(unsigned ci, unsigned j):
-        m_ci(static_cast<int>(ci)),
-        m_j(static_cast<int>(j))
-    {}
-    bool operator==(const constraint_index_and_column_struct & a) const {  return a.m_ci == m_ci && a.m_j == m_j; }
-    bool operator!=(const constraint_index_and_column_struct & a) const {  return ! (*this == a);}
-};
 
 class lar_solver : public column_namer {
     //////////////////// fields //////////////////////////
     lp_settings m_settings;
-    stacked_value<lp_status> m_status = OPTIMAL;
+    stacked_value<lp_status> m_status;
-	std::unordered_map<unsigned, var_index> m_ext_vars_to_columns;
+    stacked_value<simplex_strategy_enum> m_simplex_strategy;
+    std::unordered_map<unsigned, var_index> m_ext_vars_to_columns;
     vector<unsigned> m_columns_to_ext_vars_or_term_indices;
     stacked_vector<ul_pair> m_vars_to_ul_pairs;
     vector<lar_base_constraint*> m_constraints;
     stacked_value<unsigned> m_constraint_count;
-    indexed_vector<mpq> m_incoming_buffer;
     // the set of column indices j such that bounds have changed for j
     int_set m_columns_with_changed_bound;
     int_set m_rows_with_changed_bounds;
-	int_set m_basic_columns_with_changed_cost;
+    int_set m_basic_columns_with_changed_cost;
-    stacked_value<int> m_infeasible_column_index = -1; // such can be found at the initialization step
+    stacked_value<int> m_infeasible_column_index; // such can be found at the initialization step
     stacked_value<unsigned> m_term_count;
-public: // debug remove later
     vector<lar_term*> m_terms;
-private:
     vector<lar_term*> m_orig_terms;
-    const var_index m_terms_start_index = 1000000;
+    const var_index m_terms_start_index;
     indexed_vector<mpq> m_column_buffer;    
-    std::function<column_type (unsigned)> m_column_type_function = [this] (unsigned j) {return m_mpq_lar_core_solver.m_column_types()[j];};
+    std::function<column_type (unsigned)> m_column_type_function;    
-    
 public:
     lar_core_solver m_mpq_lar_core_solver;
-	unsigned constraint_count() const {
+    unsigned constraint_count() const {
-		return m_constraints.size();
+        return m_constraints.size();
-	}
+    }
-	const lar_base_constraint& get_constraint(unsigned ci) const {
+    const lar_base_constraint& get_constraint(unsigned ci) const {
-		return *(m_constraints[ci]);
+        return *(m_constraints[ci]);
-	}
+    }
 
     ////////////////// methods ////////////////////////////////
     static_matrix<mpq, numeric_pair<mpq>> & A_r() { return m_mpq_lar_core_solver.m_r_A;}
@@ -128,10 +80,12 @@ public:
     }
 
 
-    lar_solver() : m_mpq_lar_core_solver(
+    lar_solver() : m_status(OPTIMAL),
-                                         m_settings,
+                   m_infeasible_column_index(-1),
-                                         *this
+                   m_terms_start_index(1000000),
-                                         ) {}
+                   m_column_type_function ([this] (unsigned j) {return m_mpq_lar_core_solver.m_column_types()[j];}),
+                   m_mpq_lar_core_solver(m_settings, *this)
+    {}
     
     void set_propagate_bounds_on_pivoted_rows_mode(bool v) {
         m_mpq_lar_core_solver.m_r_solver.m_pivoted_rows = v? (& m_rows_with_changed_bounds) : nullptr;
@@ -146,24 +100,7 @@ public:
             delete t;
     }
 
-    var_index add_var(unsigned ext_j) {
+#include "util/lp/init_lar_solver.h"
-        var_index i;
-        lean_assert (ext_j < m_terms_start_index); 
-
-		if (ext_j >= m_terms_start_index)
-			throw 0; // todo : what is the right was to exit?
-            
-        if (try_get_val(m_ext_vars_to_columns, ext_j, i)) {
-            return i;
-        }
-        lean_assert(m_vars_to_ul_pairs.size() == A_r().column_count());
-        i = A_r().column_count();
-        m_vars_to_ul_pairs.push_back (ul_pair(static_cast<unsigned>(-1)));
-        register_new_ext_var_index(ext_j);
-        add_non_basic_var_to_core_fields();
-        lean_assert(sizes_are_correct());
-        return i;
-    }
     
     numeric_pair<mpq> const& get_value(var_index vi) const { return m_mpq_lar_core_solver.m_r_x[vi]; }
 
@@ -177,98 +114,16 @@ public:
     }
 
 
-	bool need_to_presolve_with_doubles() const { return m_mpq_lar_core_solver.need_to_presolve_with_double_solver(); }
+    bool use_lu() const { return m_settings.simplex_strategy() == simplex_strategy_enum::lu; }
-
-    void add_row_from_term_no_constraint(const lar_term * term) {
-        // j will be a new variable 
-        unsigned j = A_r().column_count();
-        ul_pair ul(j);
-        m_vars_to_ul_pairs.push_back(ul);
-        add_basic_var_to_core_fields();
-        if (use_tableau()) {
-            auto it = iterator_on_term_with_basis_var(*term, j);
-            A_r().fill_last_row_with_pivoting(it,
-                                              m_mpq_lar_core_solver.m_r_solver.m_basis_heading);
-            m_mpq_lar_core_solver.m_r_solver.m_b.resize(A_r().column_count(), zero_of_type<mpq>());
-        } else {
-            fill_last_row_of_A_r(A_r(), term);
-        }
-        m_mpq_lar_core_solver.m_r_x[j] = get_basic_var_value_from_row_directly(A_r().row_count() - 1);
-		if (need_to_presolve_with_doubles())
-			fill_last_row_of_A_d(A_d(), term);
-    }
-    
-    void add_constraint_from_term_and_create_new_column_row(unsigned term_j, const lar_term* term,
-                                                            lconstraint_kind kind, const mpq & right_side) {
-
-        lean_assert(A_r().column_count() == m_mpq_lar_core_solver.m_r_solver.m_costs.size());
-        // j will be a new variable 
-        unsigned j = A_r().column_count();
-        ul_pair ul(j);
-        m_vars_to_ul_pairs.push_back(ul);
-        add_basic_var_to_core_fields();
-        if (!m_settings.use_tableau()) {
-            fill_last_row_of_A_r(A_r(), term);
-        }
-        else {
-            auto it = iterator_on_term_with_basis_var(*term, j);
-            A_r().fill_last_row_with_pivoting(it,
-                                              m_mpq_lar_core_solver.m_r_solver.m_basis_heading);
-            m_mpq_lar_core_solver.m_r_solver.m_b.resize(A_r().column_count(), zero_of_type<mpq>());
-        }
-        m_mpq_lar_core_solver.m_r_x[A_r().column_count() - 1] = get_basic_var_value_from_row_directly(A_r().row_count() - 1);
-        fill_last_row_of_A_d(A_d(), term);
-        register_new_ext_var_index(term_j);
-
-        // m_constraints.size() is the index of the constrained that is about to be added
-        update_column_type_and_bound(j, kind, right_side - term->m_v, m_constraints.size());
-        m_constraints.push_back(new lar_term_constraint(term, kind, right_side));
-        lean_assert(A_r().column_count() == m_mpq_lar_core_solver.m_r_solver.m_costs.size());
-    }
-
-	void add_var_bound_on_constraint_for_term(var_index j, lconstraint_kind kind, const mpq & right_side, constraint_index ci) {
-		lean_assert(is_term(j));
-		unsigned adjusted_term_index = adjust_term_index(j);
-		unsigned term_j;
-		if (try_get_val(m_ext_vars_to_columns, j, term_j)) {
-			mpq rs = right_side - m_orig_terms[adjusted_term_index]->m_v;
-			m_constraints.push_back(new lar_term_constraint(m_orig_terms[adjusted_term_index], kind, right_side));
-			update_column_type_and_bound(term_j, kind, rs, ci);
-		}
-		else {
-			add_constraint_from_term_and_create_new_column_row(j, m_orig_terms[adjusted_term_index], kind, right_side);
-		}
-	}
-
-
-    void add_row_for_term(const lar_term * term) {
-        lean_assert(sizes_are_correct());
-        add_row_from_term_no_constraint(term);
-        lean_assert(sizes_are_correct());
-    }
     
     bool sizes_are_correct() const {
-        lean_assert(!m_mpq_lar_core_solver.need_to_presolve_with_double_solver() || A_r().column_count() == A_d().column_count());
+        lean_assert(strategy_is_undecided() || !m_mpq_lar_core_solver.need_to_presolve_with_double_solver() || A_r().column_count() == A_d().column_count());
         lean_assert(A_r().column_count() == m_mpq_lar_core_solver.m_r_solver.m_column_types.size());
         lean_assert(A_r().column_count() == m_mpq_lar_core_solver.m_r_solver.m_costs.size());
         lean_assert(A_r().column_count() == m_mpq_lar_core_solver.m_r_x.size());
         return true;
     }
     
-    constraint_index add_var_bound(var_index j, lconstraint_kind kind, const mpq & right_side)  {
-        lean_assert(sizes_are_correct());
-        constraint_index ci = m_constraints.size();
-        if (!is_term(j)) { // j is a var
-            auto vc = new lar_var_constraint(j, kind, right_side);
-            m_constraints.push_back(vc);
-            update_column_type_and_bound(j, kind, right_side, ci);
-        } else {
-			add_var_bound_on_constraint_for_term(j, kind, right_side, ci);
-        }
-        lean_assert(sizes_are_correct());
-        return ci;
-    }
-
  
     void print_implied_bound(const implied_bound& be, std::ostream & out) const {
         out << "implied bound\n";
@@ -490,10 +345,10 @@ public:
         //        new linear_combination_iterator_on_vector<mpq>(m_terms[adjust_term_index(term_index)]->coeffs_as_vector());
     }
 
-	unsigned adjust_column_index_to_term_index(unsigned j) const {
+    unsigned adjust_column_index_to_term_index(unsigned j) const {
-		unsigned ext_var_or_term = m_columns_to_ext_vars_or_term_indices[j];
+        unsigned ext_var_or_term = m_columns_to_ext_vars_or_term_indices[j];
-		return ext_var_or_term < m_terms_start_index ? j : ext_var_or_term;
+        return ext_var_or_term < m_terms_start_index ? j : ext_var_or_term;
-	}
+    }
     
     void propagate_bounds_on_a_term(const lar_term& t, bound_propagator & bp, unsigned term_offset) {
         lean_assert(false); // not implemented
@@ -561,6 +416,9 @@ public:
     void set_status(lp_status s) {m_status = s;}
 
     lp_status find_feasible_solution() {
+        if (strategy_is_undecided())
+            decide_on_strategy_and_adjust_initial_state();
+
         m_mpq_lar_core_solver.m_r_solver.m_look_for_feasible_solution_only = true;
         return solve();
     }
@@ -599,7 +457,8 @@ public:
         return ret;
     }
     void push() {
-        lean_assert(sizes_are_correct());
+        m_simplex_strategy = m_settings.simplex_strategy();
+        m_simplex_strategy.push();
         m_status.push();
         m_vars_to_ul_pairs.push();
         m_infeasible_column_index.push();
@@ -608,7 +467,6 @@ public:
         m_term_count.push();
         m_constraint_count = m_constraints.size();
         m_constraint_count.push();
-        lean_assert(sizes_are_correct());
     }
 
     static void clean_large_elements_after_pop(unsigned n, int_set& set) {
@@ -627,8 +485,7 @@ public:
 
     
     void pop(unsigned k) {
-        lean_assert(sizes_are_correct());
+        int n_was = static_cast<int>(m_ext_vars_to_columns.size());
-		int n_was = static_cast<int>(m_ext_vars_to_columns.size());
 		m_status.pop(k);
 		m_infeasible_column_index.pop(k);
         unsigned n = m_vars_to_ul_pairs.peek_size(k);
@@ -645,7 +502,8 @@ public:
         unsigned m = A_r().row_count();
         clean_large_elements_after_pop(m, m_rows_with_changed_bounds);
         clean_inf_set_of_r_solver_after_pop();
-        lean_assert(!use_tableau() || m_mpq_lar_core_solver.m_r_solver.reduced_costs_are_correct_tableau());
+        lean_assert(m_settings.simplex_strategy() == simplex_strategy_enum::undecided ||
+			(!use_tableau()) || m_mpq_lar_core_solver.m_r_solver.reduced_costs_are_correct_tableau());
         
         
         lean_assert(ax_is_correct());
@@ -662,7 +520,9 @@ public:
         }
         m_terms.resize(m_term_count);
         m_orig_terms.resize(m_term_count);
-        lean_assert(sizes_are_correct());
+		m_simplex_strategy.pop(k);
+		m_settings.simplex_strategy() = m_simplex_strategy;
+		lean_assert(sizes_are_correct());
         lean_assert((!m_settings.use_tableau()) || m_mpq_lar_core_solver.m_r_solver.reduced_costs_are_correct_tableau());
     }
     
@@ -676,6 +536,9 @@ public:
 
     bool maximize_term_on_tableau(const vector<std::pair<mpq, var_index>> & term,
                                   impq &term_max) {
+        if (settings().simplex_strategy() == simplex_strategy_enum::undecided)
+            decide_on_strategy_and_adjust_initial_state();
+ 
         m_mpq_lar_core_solver.solve();
         if (m_mpq_lar_core_solver.m_r_solver.get_status() == UNBOUNDED)
             return false;
@@ -747,13 +610,13 @@ public:
     bool maximize_term_on_corrected_r_solver(const vector<std::pair<mpq, var_index>> & term,
                                              impq &term_max) {
         settings().backup_costs = false;
-        switch (settings().m_simplex_strategy) {
+        switch (settings().simplex_strategy()) {
         case simplex_strategy_enum::tableau_rows:
             prepare_costs_for_r_solver(term);
-            settings().m_simplex_strategy = simplex_strategy_enum::tableau_costs;
+            settings().simplex_strategy() = simplex_strategy_enum::tableau_costs;
             {
                 bool ret = maximize_term_on_tableau(term, term_max);
-                settings().m_simplex_strategy = simplex_strategy_enum::tableau_rows;
+                settings().simplex_strategy() = simplex_strategy_enum::tableau_rows;
                 set_costs_to_zero(term);
                 m_mpq_lar_core_solver.m_r_solver.set_status(OPTIMAL);
                 return ret;
@@ -767,7 +630,7 @@ public:
                 return ret;
             }
             
-        case simplex_strategy_enum::no_tableau:
+        case simplex_strategy_enum::lu:
             lean_assert(false); // not implemented
             return false;
         default:
@@ -786,24 +649,6 @@ public:
     
 
     
-    var_index add_term(const vector<std::pair<mpq, var_index>> & coeffs,
-                       const mpq &m_v) {
-
-        lean_assert(A_r().column_count() == m_mpq_lar_core_solver.m_r_solver.m_costs.size());
-        m_terms.push_back(new lar_term(coeffs, m_v));
-        m_orig_terms.push_back(new lar_term(coeffs, m_v));
-	    unsigned adjusted_term_index = m_terms.size() - 1;
-        if (use_tableau() && !coeffs.empty()) {
-            register_new_ext_var_index(m_terms_start_index + adjusted_term_index);
-            add_row_for_term(m_orig_terms.back());
-            if (m_settings.bound_propagation())
-                m_rows_with_changed_bounds.insert(A_r().row_count() - 1);
-        } 
-		lean_assert(m_ext_vars_to_columns.size() == A_r().column_count());
-        lean_assert(A_r().column_count() == m_mpq_lar_core_solver.m_r_solver.m_costs.size());
-        return m_terms_start_index + adjusted_term_index;
-    }
-
     const lar_term &  get_term(unsigned j) const {
         lean_assert(j >= m_terms_start_index);
         return *m_terms[j - m_terms_start_index];
@@ -818,78 +663,6 @@ public:
         A_d().pop(k);
     }
 
-    void add_new_var_to_core_fields_for_mpq(bool register_in_basis) {
-        unsigned j = A_r().column_count();
-        A_r().add_column();
-        lean_assert(m_mpq_lar_core_solver.m_r_x.size() == j);
-        //        lean_assert(m_mpq_lar_core_solver.m_r_low_bounds.size() == j && m_mpq_lar_core_solver.m_r_upper_bounds.size() == j);  // restore later
-        m_mpq_lar_core_solver.m_r_x.resize(j + 1);
-        m_mpq_lar_core_solver.m_r_low_bounds.increase_size_by_one();
-        m_mpq_lar_core_solver.m_r_upper_bounds.increase_size_by_one();
-        m_mpq_lar_core_solver.m_r_solver.m_inf_set.increase_size_by_one();
-        m_mpq_lar_core_solver.m_r_solver.m_costs.resize(j + 1);
-        m_mpq_lar_core_solver.m_r_solver.m_d.resize(j + 1);
-        lean_assert(m_mpq_lar_core_solver.m_r_heading.size() == j); // as A().column_count() on the entry to the method
-        if (register_in_basis) {
-            A_r().add_row();
-            m_mpq_lar_core_solver.m_r_heading.push_back(m_mpq_lar_core_solver.m_r_basis.size());
-            m_mpq_lar_core_solver.m_r_basis.push_back(j);
-            if (m_settings.bound_propagation())
-                m_rows_with_changed_bounds.insert(A_r().row_count() - 1);
-        } else {
-            m_mpq_lar_core_solver.m_r_heading.push_back(- static_cast<int>(m_mpq_lar_core_solver.m_r_nbasis.size()) - 1);
-            m_mpq_lar_core_solver.m_r_nbasis.push_back(j);
-        }
-    }
-
-    void add_new_var_to_core_fields_for_doubles(bool register_in_basis) {
-        unsigned j = A_d().column_count();
-        A_d().add_column();
-        lean_assert(m_mpq_lar_core_solver.m_d_x.size() == j);
-        //        lean_assert(m_mpq_lar_core_solver.m_d_low_bounds.size() == j && m_mpq_lar_core_solver.m_d_upper_bounds.size() == j);  // restore later
-        m_mpq_lar_core_solver.m_d_x.resize(j + 1 ); 
-        m_mpq_lar_core_solver.m_d_low_bounds.resize(j + 1);
-        m_mpq_lar_core_solver.m_d_upper_bounds.resize(j + 1);
-        lean_assert(m_mpq_lar_core_solver.m_d_heading.size() == j); // as A().column_count() on the entry to the method
-        if (register_in_basis) {
-            A_d().add_row();
-            m_mpq_lar_core_solver.m_d_heading.push_back(m_mpq_lar_core_solver.m_d_basis.size());
-            m_mpq_lar_core_solver.m_d_basis.push_back(j);
-        }else {
-            m_mpq_lar_core_solver.m_d_heading.push_back(- static_cast<int>(m_mpq_lar_core_solver.m_d_nbasis.size()) - 1);
-            m_mpq_lar_core_solver.m_d_nbasis.push_back(j);
-        }
-    }
-
-    void add_basic_var_to_core_fields() {
-		bool need_to_presolve_with_doubles = m_mpq_lar_core_solver.need_to_presolve_with_double_solver();
-        lean_assert(!need_to_presolve_with_doubles || A_r().column_count() == A_d().column_count());
-        m_mpq_lar_core_solver.m_column_types.push_back(column_type::free_column);
-        m_columns_with_changed_bound.increase_size_by_one();
-        m_rows_with_changed_bounds.increase_size_by_one();
-        add_new_var_to_core_fields_for_mpq(true);
-		if (need_to_presolve_with_doubles)
-			add_new_var_to_core_fields_for_doubles(true);
-    }
-
-    void add_non_basic_var_to_core_fields() {
-        lean_assert(!m_mpq_lar_core_solver.need_to_presolve_with_double_solver() || A_r().column_count() == A_d().column_count());
-        m_mpq_lar_core_solver.m_column_types.push_back(column_type::free_column);
-        m_columns_with_changed_bound.increase_size_by_one();
-        add_new_var_to_core_fields_for_mpq(false);
-		if (m_mpq_lar_core_solver.need_to_presolve_with_double_solver())
-			add_new_var_to_core_fields_for_doubles(false);
-    }
-    
-    void register_new_ext_var_index(unsigned s) {
-        lean_assert(!contains(m_ext_vars_to_columns, s));
-        unsigned j = static_cast<unsigned>(m_ext_vars_to_columns.size());
-        m_ext_vars_to_columns[s] = j;
-        lean_assert(m_columns_to_ext_vars_or_term_indices.size() == j);
-        m_columns_to_ext_vars_or_term_indices.push_back(s);
-    }
-
-    
 
     void set_upper_bound_witness(var_index j, constraint_index ci) {
         ul_pair ul = m_vars_to_ul_pairs[j];
@@ -903,312 +676,6 @@ public:
         m_vars_to_ul_pairs[j] = ul;
     }
 
-    void update_free_column_type_and_bound(var_index j, lconstraint_kind kind, const mpq & right_side, constraint_index constr_ind) {
-        mpq y_of_bound(0);
-        switch (kind) {
-        case LT:
-            y_of_bound = -1;
-        case LE:
-            m_mpq_lar_core_solver.m_column_types[j] = column_type::upper_bound;
-            lean_assert(m_mpq_lar_core_solver.m_column_types()[j] == column_type::upper_bound);
-            lean_assert(m_mpq_lar_core_solver.m_r_upper_bounds.size() > j);
-            {
-                auto up = numeric_pair<mpq>(right_side, y_of_bound);
-                m_mpq_lar_core_solver.m_r_upper_bounds[j] = up;
-            }
-            set_upper_bound_witness(j, constr_ind);
-            break;
-        case GT:
-            y_of_bound = 1;
-        case GE:
-            m_mpq_lar_core_solver.m_column_types[j] = column_type::low_bound;
-            lean_assert(m_mpq_lar_core_solver.m_r_upper_bounds.size() > j);
-            {
-                auto low = numeric_pair<mpq>(right_side, y_of_bound);
-                m_mpq_lar_core_solver.m_r_low_bounds[j] = low;
-            }
-            set_low_bound_witness(j, constr_ind);
-            break;
-        case EQ:
-            m_mpq_lar_core_solver.m_column_types[j] = column_type::fixed;
-            m_mpq_lar_core_solver.m_r_low_bounds[j] = m_mpq_lar_core_solver.m_r_upper_bounds[j] = numeric_pair<mpq>(right_side, zero_of_type<mpq>());
-            set_upper_bound_witness(j, constr_ind);
-            set_low_bound_witness(j, constr_ind);
-            break;
-
-        default:
-            lean_unreachable();
-                
-        }
-        m_columns_with_changed_bound.insert(j);
-    }
-
-    void update_upper_bound_column_type_and_bound(var_index j, lconstraint_kind kind, const mpq & right_side, constraint_index ci) {
-        lean_assert(m_mpq_lar_core_solver.m_column_types()[j] == column_type::upper_bound);
-        mpq y_of_bound(0);
-        switch (kind) {
-        case LT:
-            y_of_bound = -1;
-        case LE:
-            {
-                auto up = numeric_pair<mpq>(right_side, y_of_bound);
-                if (up < m_mpq_lar_core_solver.m_r_upper_bounds()[j]) {
-                    m_mpq_lar_core_solver.m_r_upper_bounds[j] = up;
-                    set_upper_bound_witness(j, ci);
-                    m_columns_with_changed_bound.insert(j);
-                }
-            }
-            break;
-        case GT:
-            y_of_bound = 1;
-        case GE:            
-            m_mpq_lar_core_solver.m_column_types[j] = column_type::boxed;
-            {
-                auto low = numeric_pair<mpq>(right_side, y_of_bound);
-                m_mpq_lar_core_solver.m_r_low_bounds[j] = low;
-                set_low_bound_witness(j, ci);
-                m_columns_with_changed_bound.insert(j);
-                if (low > m_mpq_lar_core_solver.m_r_upper_bounds[j]) {
-                    m_status = INFEASIBLE;
-                    m_infeasible_column_index = j;
-                } else {
-                    m_mpq_lar_core_solver.m_column_types[j] = m_mpq_lar_core_solver.m_r_low_bounds()[j] < m_mpq_lar_core_solver.m_r_upper_bounds()[j]? column_type::boxed : column_type::fixed;
-                }                     
-            }
-            break;
-        case EQ:
-            {
-                auto v = numeric_pair<mpq>(right_side, zero_of_type<mpq>());
-                if (v > m_mpq_lar_core_solver.m_r_upper_bounds[j]) {
-                    m_status = INFEASIBLE;
-                    set_low_bound_witness(j, ci);
-                    m_infeasible_column_index = j;
-                } else {
-                    m_mpq_lar_core_solver.m_r_low_bounds[j] = m_mpq_lar_core_solver.m_r_upper_bounds[j] = v;
-                    m_columns_with_changed_bound.insert(j);
-                    set_low_bound_witness(j, ci);
-                    set_upper_bound_witness(j, ci);
-                    m_mpq_lar_core_solver.m_column_types[j] = column_type::fixed;
-                }
-                break;
-            }
-            break;
-
-        default:
-            lean_unreachable();
-                
-        }
-    }
-    
-    void update_boxed_column_type_and_bound(var_index j, lconstraint_kind kind, const mpq & right_side, constraint_index ci) {
-        lean_assert(m_status == INFEASIBLE || (m_mpq_lar_core_solver.m_column_types()[j] == column_type::boxed && m_mpq_lar_core_solver.m_r_low_bounds()[j] < m_mpq_lar_core_solver.m_r_upper_bounds()[j]));
-        mpq y_of_bound(0);
-        switch (kind) {
-        case LT:
-            y_of_bound = -1;
-        case LE:
-            {
-                auto up = numeric_pair<mpq>(right_side, y_of_bound);
-                if (up < m_mpq_lar_core_solver.m_r_upper_bounds[j]) {
-                    m_mpq_lar_core_solver.m_r_upper_bounds[j] = up;
-                    set_upper_bound_witness(j, ci);
-                    m_columns_with_changed_bound.insert(j);
-                }
-
-                if (up < m_mpq_lar_core_solver.m_r_low_bounds[j]) {
-                    m_status = INFEASIBLE;
-                    lean_assert(false);
-                    m_infeasible_column_index = j;
-                } else {
-                    if (m_mpq_lar_core_solver.m_r_low_bounds()[j] == m_mpq_lar_core_solver.m_r_upper_bounds()[j])
-                        m_mpq_lar_core_solver.m_column_types[j] = column_type::fixed;
-                }                    
-            }
-            break;
-        case GT:
-            y_of_bound = 1;
-        case GE:            
-            {
-                auto low = numeric_pair<mpq>(right_side, y_of_bound);
-                if (low > m_mpq_lar_core_solver.m_r_low_bounds[j]) {
-                    m_mpq_lar_core_solver.m_r_low_bounds[j] = low;
-                    m_columns_with_changed_bound.insert(j);
-                    set_low_bound_witness(j, ci);
-                }
-                if (low > m_mpq_lar_core_solver.m_r_upper_bounds[j]) {
-                    m_status = INFEASIBLE;
-                    m_infeasible_column_index = j;
-                } else if ( low == m_mpq_lar_core_solver.m_r_upper_bounds[j]) {
-                    m_mpq_lar_core_solver.m_column_types[j] = column_type::fixed;
-                }
-            }
-            break;
-        case EQ:
-            {
-                auto v = numeric_pair<mpq>(right_side, zero_of_type<mpq>());
-                if (v < m_mpq_lar_core_solver.m_r_low_bounds[j]) {
-                    m_status = INFEASIBLE;
-                    m_infeasible_column_index = j;
-                    set_upper_bound_witness(j, ci);                    
-                } else if (v > m_mpq_lar_core_solver.m_r_upper_bounds[j]) {
-                    m_status = INFEASIBLE;
-                    m_infeasible_column_index = j;
-                    set_low_bound_witness(j, ci);                    
-                } else {
-                    m_mpq_lar_core_solver.m_r_low_bounds[j] = m_mpq_lar_core_solver.m_r_upper_bounds[j] = v;
-                    set_low_bound_witness(j, ci);
-                    set_upper_bound_witness(j, ci);
-                    m_mpq_lar_core_solver.m_column_types[j] = column_type::fixed;
-                    m_columns_with_changed_bound.insert(j);
-                }
-                
-                break;
-            }
-
-        default:
-            lean_unreachable();
-                
-        }
-    }
-    void update_low_bound_column_type_and_bound(var_index j, lconstraint_kind kind, const mpq & right_side, constraint_index ci) {
-        lean_assert(m_mpq_lar_core_solver.m_column_types()[j] == column_type::low_bound);
-        mpq y_of_bound(0);
-        switch (kind) {
-        case LT:
-            y_of_bound = -1;
-        case LE:
-            {
-                auto up = numeric_pair<mpq>(right_side, y_of_bound);
-                m_mpq_lar_core_solver.m_r_upper_bounds[j] = up;
-                set_upper_bound_witness(j, ci);
-                m_columns_with_changed_bound.insert(j);
-
-                if (up < m_mpq_lar_core_solver.m_r_low_bounds[j]) {
-                    m_status = INFEASIBLE;
-                    m_infeasible_column_index = j;
-                } else {
-                    m_mpq_lar_core_solver.m_column_types[j] = m_mpq_lar_core_solver.m_r_low_bounds()[j] < m_mpq_lar_core_solver.m_r_upper_bounds()[j]? column_type::boxed : column_type::fixed;
-                }                    
-            }
-            break;
-        case GT:
-            y_of_bound = 1;
-        case GE:            
-            {
-                auto low = numeric_pair<mpq>(right_side, y_of_bound);
-                if (low > m_mpq_lar_core_solver.m_r_low_bounds[j]) {
-                    m_mpq_lar_core_solver.m_r_low_bounds[j] = low;
-                    m_columns_with_changed_bound.insert(j);
-                    set_low_bound_witness(j, ci);
-                }
-            }
-            break;
-        case EQ:
-            {
-                auto v = numeric_pair<mpq>(right_side, zero_of_type<mpq>());
-                if (v < m_mpq_lar_core_solver.m_r_low_bounds[j]) {
-                    m_status = INFEASIBLE;
-                    m_infeasible_column_index = j;
-                    set_upper_bound_witness(j, ci);                    
-                } else {
-                    m_mpq_lar_core_solver.m_r_low_bounds[j] = m_mpq_lar_core_solver.m_r_upper_bounds[j] = v;
-                    set_low_bound_witness(j, ci);
-                    set_upper_bound_witness(j, ci);
-                    m_mpq_lar_core_solver.m_column_types[j] = column_type::fixed;
-                }
-                m_columns_with_changed_bound.insert(j);
-                break;
-            }
-
-        default:
-            lean_unreachable();
-                
-        }
-    }
-
-    void update_fixed_column_type_and_bound(var_index j, lconstraint_kind kind, const mpq & right_side, constraint_index ci) {
-        lean_assert(m_status == INFEASIBLE || (m_mpq_lar_core_solver.m_column_types()[j] == column_type::fixed && m_mpq_lar_core_solver.m_r_low_bounds()[j] == m_mpq_lar_core_solver.m_r_upper_bounds()[j]));
-        lean_assert(m_status == INFEASIBLE || (m_mpq_lar_core_solver.m_r_low_bounds()[j].y.is_zero() && m_mpq_lar_core_solver.m_r_upper_bounds()[j].y.is_zero()));
-        auto v = numeric_pair<mpq>(right_side, mpq(0));
-        
-        mpq y_of_bound(0);
-        switch (kind) {
-        case LT:
-            if (v <= m_mpq_lar_core_solver.m_r_low_bounds[j]) {
-                m_status = INFEASIBLE;
-                m_infeasible_column_index = j;
-                set_upper_bound_witness(j, ci);
-            }                   
-            break;
-        case LE:
-            {
-                if (v < m_mpq_lar_core_solver.m_r_low_bounds[j]) {
-                    m_status = INFEASIBLE;
-                    m_infeasible_column_index = j;
-                    set_upper_bound_witness(j, ci);
-                }                   
-            }
-            break;
-        case GT:
-            {
-                if (v >= m_mpq_lar_core_solver.m_r_upper_bounds[j]) {
-                    m_status = INFEASIBLE;
-                    m_infeasible_column_index =j;
-                    set_low_bound_witness(j, ci);
-                }
-            }
-            break;
-        case GE:            
-            {
-                if (v > m_mpq_lar_core_solver.m_r_upper_bounds[j]) {
-                    m_status = INFEASIBLE;
-                    m_infeasible_column_index = j;
-                    set_low_bound_witness(j, ci);
-                }
-            }
-            break;
-        case EQ:
-            {
-                if (v < m_mpq_lar_core_solver.m_r_low_bounds[j]) {
-                    m_status = INFEASIBLE;
-                    m_infeasible_column_index = j;
-                    set_upper_bound_witness(j, ci);                    
-                } else if (v > m_mpq_lar_core_solver.m_r_upper_bounds[j]) {
-                    m_status = INFEASIBLE;
-                    m_infeasible_column_index = j;
-                    set_low_bound_witness(j, ci);                    
-                } 
-                break;
-            }
-
-        default:
-            lean_unreachable();
-                
-        }
-    }
-    
-    void update_column_type_and_bound(var_index j, lconstraint_kind kind, const mpq & right_side, constraint_index constr_index) {
-        switch(m_mpq_lar_core_solver.m_column_types[j]) {
-        case column_type::free_column:
-            update_free_column_type_and_bound(j, kind, right_side, constr_index);
-            break;
-        case column_type::boxed:
-            update_boxed_column_type_and_bound(j, kind, right_side, constr_index);
-            break;
-        case column_type::low_bound:
-            update_low_bound_column_type_and_bound(j, kind, right_side, constr_index);
-            break;
-        case column_type::upper_bound:
-            update_upper_bound_column_type_and_bound(j, kind, right_side, constr_index);
-            break;
-        case column_type::fixed:
-            update_fixed_column_type_and_bound(j, kind, right_side, constr_index);
-            break;
-        default:
-            lean_assert(false); // cannot be here
-        }
-    }
-
 
     void substitute_terms(const mpq & mult,
                           const vector<std::pair<mpq, var_index>>& left_side_with_terms,
@@ -1247,8 +714,9 @@ public:
 
     bool use_tableau() const { return m_settings.use_tableau(); }
 
-    bool use_tableau_costs() const { return m_settings.simplex_strategy() == simplex_strategy_enum::tableau_costs; }
+    bool use_tableau_costs() const {
-
+        return m_settings.simplex_strategy() == simplex_strategy_enum::tableau_costs;
+    }
     
     void detect_rows_of_column_with_bound_change(unsigned j) {
         if (m_mpq_lar_core_solver.m_r_heading[j] >= 0) { // it is a basic column
@@ -1485,23 +953,6 @@ public:
         unsigned basis_j = A.column_count() - 1;
         A.set(last_row, basis_j, mpq(1));
     }
-    // this fills the last row of A_d and sets the basis column: -1 in the last column of the row
-    void fill_last_row_of_A_d(static_matrix<double, double> & A, const lar_term* ls) {
-        lean_assert(A.row_count() > 0);
-        lean_assert(A.column_count() > 0);
-        unsigned last_row = A.row_count() - 1;
-        lean_assert(A.m_rows[last_row].empty());
-
-        for (auto & t : ls->m_coeffs) {
-            lean_assert(!is_zero(t.second));
-            var_index j = t.first;
-            A.set(last_row, j, - t.second.get_double());
-        }
-
-        unsigned basis_j = A.column_count() - 1;
-        A.set(last_row, basis_j, - 1 );
-    }
-
 
     template <typename U, typename V>
     void create_matrix_A(static_matrix<U, V> & matr) {
@@ -1518,8 +969,8 @@ public:
 
     template <typename U, typename V>
     void copy_from_mpq_matrix(static_matrix<U, V> & matr) {
-        lean_assert(matr.row_count() == A_r().row_count());
+		matr.m_rows.resize(A_r().row_count());
-        lean_assert(matr.column_count() == A_r().column_count());
+		matr.m_columns.resize(A_r().column_count());
         for (unsigned i = 0; i < matr.row_count(); i++) {
             for (auto & it : A_r().m_rows[i]) {
                 matr.set(i, it.m_j,  convert_struct<U, mpq>::convert(it.get_val()));
@@ -1691,10 +1142,8 @@ public:
 
     bool explanation_is_correct(const vector<std::pair<mpq, unsigned>>& explanation) const {
 #ifdef LEAN_DEBUG
-#if 0
-        // disabled as 'kind' is not assigned
         lconstraint_kind kind;
-        the_relations_are_of_same_type(explanation, kind);
+        lean_assert(the_relations_are_of_same_type(explanation, kind));
         lean_assert(the_left_sides_sum_to_zero(explanation));
         mpq rs = sum_of_right_sides_of_explanation(explanation);
         switch (kind) {
@@ -1712,7 +1161,6 @@ public:
             lean_assert(false);
             return false;
         }
-#endif
 #endif
         return true;
     }
@@ -1737,58 +1185,6 @@ public:
         return ret;
     }
 
-    // template <typename U, typename V>
-    // void prepare_core_solver_fields_with_signature(static_matrix<U, V> & A, vector<V> & x,
-    //                                                            vector<V> & low_bound,
-    //                                                            vector<V> & upper_bound, const lar_solution_signature & signature) {
-    //     create_matrix_A_r(A);
-    //     fill_bounds_for_core_solver(low_bound, upper_bound);
-    //     if (m_status == INFEASIBLE) {
-    //         lean_assert(false); // not implemented
-    //     }
-
-    //     resize_and_init_x_with_signature(x, low_bound, upper_bound, signature);
-    //     lean_assert(A.column_count() == x.size());
-    // }
-
-    // void find_solution_signature_with_doubles(lar_solution_signature & signature) {
-    //     static_matrix<double, double> A;
-    //     vector<double> x, low_bounds, upper_bounds;
-    //     lean_assert(false); // not implemented
-    //     // prepare_core_solver_fields<double, double>(A, x, low_bounds, upper_bounds);
-    //     vector<double> column_scale_vector;
-    //     vector<double> right_side_vector(A.row_count(), 0);
-
-    //     scaler<double, double > scaler(right_side_vector,
-    //         A,
-    //         m_settings.scaling_minimum,
-    //         m_settings.scaling_maximum,
-    //         column_scale_vector,
-    //         m_settings);
-    //     if (!scaler.scale()) {
-    //         // the scale did not succeed, unscaling
-    //         A.clear();
-    //         create_matrix_A_r(A);
-    //         for (auto & s : column_scale_vector)
-    //             s = 1.0;
-    //     }
-    //     vector<double> costs(A.column_count());
-    //     auto core_solver = lp_primal_core_solver<double, double>(A,
-    //                                                              right_side_vector,
-    //                                                              x,
-    //                                                              m_mpq_lar_core_solver.m_basis,
-    //                                                              m_mpq_lar_core_solver.m_nbasis,
-    //                                                              m_mpq_lar_core_solver.m_heading,
-    //                                                              costs,
-    //                                                              m_mpq_lar_core_solver.m_column_types(),
-    //                                                              low_bounds,
-    //                                                              upper_bounds,
-    //                                                              m_settings,
-    //                                                              *this);
-    //     core_solver.find_feasible_solution();
-    //     extract_signature_from_lp_core_solver(core_solver, signature);
-    // }
-
     bool has_lower_bound(var_index var, constraint_index& ci, mpq& value, bool& is_strict) {
 
         if (var >= m_vars_to_ul_pairs.size()) {
@@ -1865,12 +1261,29 @@ public:
 
 
     void get_model(std::unordered_map<var_index, mpq> & variable_values) const {
+        mpq delta = mpq(1, 2); // start from 0.5 to have less clashes
         lean_assert(m_status == OPTIMAL);
-        mpq delta = m_mpq_lar_core_solver.find_delta_for_strict_bounds();
+        unsigned i;
-        for (unsigned i = 0; i < m_mpq_lar_core_solver.m_r_x.size(); i++ ) {
+        do {
-            const numeric_pair<mpq> & rp = m_mpq_lar_core_solver.m_r_x[i];
+            
-            variable_values[i] = rp.x + delta * rp.y;
+            // different pairs have to produce different singleton values
-        }
+            std::unordered_set<impq> set_of_different_pairs; 
+            std::unordered_set<mpq> set_of_different_singles;
+            delta = m_mpq_lar_core_solver.find_delta_for_strict_bounds(delta);
+            for (i = 0; i < m_mpq_lar_core_solver.m_r_x.size(); i++ ) {
+                const numeric_pair<mpq> & rp = m_mpq_lar_core_solver.m_r_x[i];
+                set_of_different_pairs.insert(rp);
+                mpq x = rp.x + delta * rp.y;
+                set_of_different_singles.insert(x);
+                if (set_of_different_pairs.size()
+                    != set_of_different_singles.size()) {
+                    delta /= mpq(2);
+                    break;
+                }
+                    
+                variable_values[i] = x;
+            }
+        } while (i != m_mpq_lar_core_solver.m_r_x.size());
     }
 
 

src/util/lp/linear_combination_iterator.h

@@ -16,7 +16,7 @@ struct linear_combination_iterator {
 template <typename T>
 struct linear_combination_iterator_on_vector : linear_combination_iterator<T> {
     vector<std::pair<T, unsigned>> & m_vector;
-    int m_offset = 0;
+    int m_offset;
     bool next(T & a, unsigned & i) {
         if(m_offset >= m_vector.size())
             return false;
@@ -40,7 +40,10 @@ struct linear_combination_iterator_on_vector : linear_combination_iterator<T> {
     linear_combination_iterator<T> * clone() {
         return new linear_combination_iterator_on_vector(m_vector);
     }
-    linear_combination_iterator_on_vector(vector<std::pair<T, unsigned>> & vec): m_vector(vec) {}
+    linear_combination_iterator_on_vector(vector<std::pair<T, unsigned>> & vec):
+        m_vector(vec),
+        m_offset(0)
+    {}
     unsigned size() const { return m_vector.size(); }
 };
 

src/util/lp/lp_core_solver_base.h

@@ -17,7 +17,8 @@ namespace lean {
 
 template <typename T, typename X> // X represents the type of the x variable and the bounds
 class lp_core_solver_base {    
-    unsigned m_total_iterations = 0;
+    unsigned m_total_iterations;
+    unsigned m_iters_with_no_cost_growing;
     unsigned inc_total_iterations() { ++m_settings.st().m_total_iterations; return m_total_iterations++; }
 private:
     lp_status m_status;
@@ -25,40 +26,39 @@ public:
     bool current_x_is_feasible() const { return m_inf_set.size() == 0; }
     bool current_x_is_infeasible() const { return m_inf_set.size() != 0; }
     int_set m_inf_set;
-    bool m_using_infeas_costs = false;
+    bool m_using_infeas_costs;
 
 
-    vector<unsigned> m_columns_nz; // m_columns_nz[i] keeps an approximate value of non zeroes the i-th column
+    vector<unsigned>      m_columns_nz; // m_columns_nz[i] keeps an approximate value of non zeroes the i-th column
-    vector<unsigned> m_rows_nz; // m_rows_nz[i] keeps an approximate value of non zeroes in the i-th row
+    vector<unsigned>      m_rows_nz; // m_rows_nz[i] keeps an approximate value of non zeroes in the i-th row
-    indexed_vector<T> m_pivot_row_of_B_1;  // the pivot row of the reverse of B
+    indexed_vector<T>     m_pivot_row_of_B_1;  // the pivot row of the reverse of B
-    indexed_vector<T> m_pivot_row; // this is the real pivot row of the simplex tableu
+    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> & m_b; // the right side
+    vector<X> &           m_b; // the right side
-    vector<unsigned> & m_basis;
+    vector<unsigned> &    m_basis;
-    vector<unsigned>& m_nbasis;
+    vector<unsigned>&     m_nbasis;
-    vector<int>& m_basis_heading;
+    vector<int>&          m_basis_heading;
-    vector<X> & m_x; // a feasible solution, the fist time set in the constructor
+    vector<X> &           m_x; // a feasible solution, the fist time set in the constructor
-    vector<T> & m_costs;
+    vector<T> &           m_costs;
-    lp_settings & m_settings;
+    lp_settings &         m_settings;
-    vector<T> m_y; // the buffer for yB = cb
+    vector<T>             m_y; // the buffer for yB = cb
     // a device that is able to solve Bx=c, xB=d, and change the basis
-    lu<T, X> * m_factorization = nullptr;
+    lu<T, X> *            m_factorization;
-    const column_namer & m_column_names;
+    const column_namer &  m_column_names;
-    indexed_vector<T> m_w; // the vector featuring in 24.3 of the Chvatal book
+    indexed_vector<T>     m_w; // the vector featuring in 24.3 of the Chvatal book
-    vector<T> m_d; // the vector of reduced costs
+    vector<T>             m_d; // the vector of reduced costs
-    indexed_vector<T> m_ed; // the solution of B*m_ed = a
+    indexed_vector<T>     m_ed; // the solution of B*m_ed = a
-    unsigned m_iters_with_no_cost_growing = 0;
     const vector<column_type> & m_column_types;
-    const vector<X> & m_low_bounds;
+    const vector<X> &     m_low_bounds;
-    const vector<X> & m_upper_bounds;
+    const vector<X> &     m_upper_bounds;
-    vector<T> m_column_norms; // the approximate squares of column norms that help choosing a profitable column
+    vector<T>             m_column_norms; // the approximate squares of column norms that help choosing a profitable column
-    vector<X> m_copy_of_xB;
+    vector<X>             m_copy_of_xB;
-    unsigned m_basis_sort_counter = 0;
+    unsigned              m_basis_sort_counter;
-    vector<T> m_steepest_edge_coefficients;
+    vector<T>             m_steepest_edge_coefficients;
-    vector<unsigned> m_trace_of_basis_change_vector; // the even positions are entering, the odd positions are leaving
+    vector<unsigned>      m_trace_of_basis_change_vector; // the even positions are entering, the odd positions are leaving
-    bool m_tracing_basis_changes = false;
+    bool                  m_tracing_basis_changes;
-    int_set* m_pivoted_rows = nullptr;
+    int_set*              m_pivoted_rows;
-    bool m_look_for_feasible_solution_only = false;
+    bool                  m_look_for_feasible_solution_only;
     void start_tracing_basis_changes() {
         m_trace_of_basis_change_vector.resize(0);
         m_tracing_basis_changes = true;
@@ -348,7 +348,7 @@ public:
             if (x_is_at_bound(j))
                 break; // we should preserve x if possible
             // snap randomly
-            if (my_random() % 2 == 1) 
+            if (m_settings.random_next() % 2 == 1) 
                 m_x[j] = m_low_bounds[j];
             else
                 m_x[j] = m_upper_bounds[j];
@@ -678,6 +678,13 @@ public:
                 lean_assert(is_zero(this->m_costs[j]));
         }
         return true;
-}
+    }
+    unsigned & iters_with_no_cost_growing() {
+        return m_iters_with_no_cost_growing;
+    }
+
+    const unsigned & iters_with_no_cost_growing() const {
+        return m_iters_with_no_cost_growing;
+    }
 };
 }

src/util/lp/lp_core_solver_base.hpp

@@ -22,8 +22,11 @@ lp_core_solver_base(static_matrix<T, X> & A,
                     const vector<column_type> & column_types,
                     const vector<X> & low_bound_values,
                     const vector<X> & upper_bound_values):
+    m_total_iterations(0),
+    m_iters_with_no_cost_growing(0),
     m_status(FEASIBLE),
     m_inf_set(A.column_count()),
+    m_using_infeas_costs(false),
     m_pivot_row_of_B_1(A.row_count()),
     m_pivot_row(A.column_count()),
     m_A(A),
@@ -45,7 +48,11 @@ lp_core_solver_base(static_matrix<T, X> & A,
     m_upper_bounds(upper_bound_values),
     m_column_norms(m_n()),
     m_copy_of_xB(m_m()),
-    m_steepest_edge_coefficients(A.column_count()) {
+    m_basis_sort_counter(0),
+    m_steepest_edge_coefficients(A.column_count()),
+    m_tracing_basis_changes(false),
+    m_pivoted_rows(nullptr),
+    m_look_for_feasible_solution_only(false) {
     lean_assert(bounds_for_boxed_are_set_correctly());    
     init();
     init_basis_heading_and_non_basic_columns_vector();
@@ -58,9 +65,8 @@ allocate_basis_heading() { // the rest of initilization will be handled by the f
 }
 template <typename T, typename X> void lp_core_solver_base<T, X>::
 init() {    
-    my_random_init(m_settings.random_seed);
     allocate_basis_heading();
-    if (!use_tableau())
+    if (m_settings.use_lu())
         init_factorization(m_factorization, m_A, m_basis, m_settings);
 }
 
@@ -527,13 +533,19 @@ update_basis_and_x(int entering, int leaving, X const & tt) {
             init_factorization(m_factorization, m_A, m_basis, m_settings);
             if (!find_x_by_solving()) {
                 restore_x(entering, tt);
-                lean_assert(!A_mult_x_is_off());
+                if(A_mult_x_is_off()) {
+                    m_status = FLOATING_POINT_ERROR;
+                    m_iters_with_no_cost_growing++;
+                    return false;
+                }
+                    
                 init_factorization(m_factorization, m_A, m_basis, m_settings);
                 m_iters_with_no_cost_growing++;
                 if (m_factorization->get_status() != LU_status::OK) {
                     std::stringstream s;
-                    s << "failing refactor on off_result for entering = " << entering << ", leaving = " << leaving << " total_iterations = " << total_iterations();
+                    //                    s << "failing refactor on off_result for entering = " << entering << ", leaving = " << leaving << " total_iterations = " << total_iterations();
-                    throw_exception(s.str());
+                    m_status = FLOATING_POINT_ERROR;
+                    return false;
                 }
                 return false;
             }

src/util/lp/lp_dual_core_solver.hpp

@@ -429,7 +429,7 @@ template <typename T, typename X> bool lp_dual_core_solver<T, X>::basis_change_a
     if (snap_runaway_nonbasic_column(m_p)) {
         if (!this->find_x_by_solving()) {
             revert_to_previous_basis();
-            this->m_iters_with_no_cost_growing++;
+            this->iters_with_no_cost_growing()++;
             return false;
         }
     }
@@ -437,7 +437,7 @@ template <typename T, typename X> bool lp_dual_core_solver<T, X>::basis_change_a
     if (!problem_is_dual_feasible()) {
         // todo : shift the costs!!!!
         revert_to_previous_basis();
-        this->m_iters_with_no_cost_growing++;
+        this->iters_with_no_cost_growing()++;
         return false;
     }
 
@@ -537,7 +537,7 @@ template <typename T, typename X> unsigned lp_dual_core_solver<T, X>::get_number
     if (this->m_m() > 300) {
         s = (unsigned)((s / 100.0) * this->m_settings.percent_of_entering_to_check);
     }
-    return my_random() % s + 1;
+    return this->m_settings.random_next() % s + 1;
 }
 
 template <typename T, typename X> bool lp_dual_core_solver<T, X>::delta_keeps_the_sign(int initial_delta_sign, const T & delta) {
@@ -715,7 +715,7 @@ template <typename T, typename X> void lp_dual_core_solver<T, X>::update_xb_afte
 
 template <typename T, typename X> void lp_dual_core_solver<T, X>::one_iteration() {
     unsigned number_of_rows_to_try = get_number_of_rows_to_try_for_leaving();
-    unsigned offset_in_rows = my_random() % this->m_m();
+    unsigned offset_in_rows = this->m_settings.random_next() % this->m_m();
     if (this->get_status() == TENTATIVE_DUAL_UNBOUNDED) {
         number_of_rows_to_try = this->m_m();
     } else {
@@ -730,14 +730,14 @@ template <typename T, typename X> void lp_dual_core_solver<T, X>::solve() { // s
     lean_assert(problem_is_dual_feasible());
     lean_assert(this->basis_heading_is_correct());
     this->set_total_iterations(0);
-    this->m_iters_with_no_cost_growing = 0;
+    this->iters_with_no_cost_growing() = 0;
     do {
         if (this->print_statistics_with_iterations_and_nonzeroes_and_cost_and_check_that_the_time_is_over("", *this->m_settings.get_message_ostream())){
             return;
         }
         one_iteration();
     } while (this->get_status() != FLOATING_POINT_ERROR && this->get_status() != DUAL_UNBOUNDED && this->get_status() != OPTIMAL &&
-             this->m_iters_with_no_cost_growing <= this->m_settings.max_number_of_iterations_with_no_improvements
+             this->iters_with_no_cost_growing() <= this->m_settings.max_number_of_iterations_with_no_improvements
              && this->total_iterations() <= this->m_settings.max_total_number_of_iterations);
 }
 }

src/util/lp/lp_dual_simplex.h

@@ -11,7 +11,7 @@ namespace lean {
 
 template <typename T, typename X>
 class lp_dual_simplex: public lp_solver<T, X> {
-    lp_dual_core_solver<T, X> * m_core_solver = nullptr;
+    lp_dual_core_solver<T, X> * m_core_solver;
     vector<T> m_b_copy;
     vector<T> m_low_bounds; // We don't have a convention here that all low bounds are zeros. At least it does not hold for the first stage solver
     vector<column_type> m_column_types_of_core_solver;
@@ -24,6 +24,7 @@ public:
         }
     }
 
+    lp_dual_simplex() : m_core_solver(nullptr) {}
 
 
     void decide_on_status_after_stage1();

src/util/lp/lp_params.pyg

@@ -3,7 +3,6 @@ def_module_params('lp',
                   params=(
                   ('rep_freq', UINT, 0, 'the report frequency, in how many iterations print the cost and other info '),
                    ('min', BOOL, False, 'minimize cost'),
-                   ('presolve_with_dbl', BOOL, True, 'presolve with double'),
                    ('print_stats', BOOL, False, 'print statistic'),
                    ('simplex_strategy', UINT, 0, 'simplex strategy for the solver'),
                    ('bprop_on_pivoted_rows', BOOL, True, 'propagate bounds on rows changed by the pivot operation')

src/util/lp/lp_primal_core_solver.h

@@ -32,21 +32,21 @@ class lp_primal_core_solver:public lp_core_solver_base<T, X> {
 public:
     // m_sign_of_entering is set to 1 if the entering variable needs
     // to grow and is set to -1  otherwise
-    unsigned m_column_norm_update_counter;
+    unsigned       m_column_norm_update_counter;
-    T m_enter_price_eps;
+    T              m_enter_price_eps;
-    int m_sign_of_entering_delta;
+    int            m_sign_of_entering_delta;
     vector<breakpoint<X>> m_breakpoints;
     binary_heap_priority_queue<X> m_breakpoint_indices_queue;
     indexed_vector<T> m_beta; // see Swietanowski working vector beta for column norms
-    T m_epsilon_of_reduced_cost = T(1)/T(10000000);
+    T                 m_epsilon_of_reduced_cost;
-    vector<T> m_costs_backup;
+    vector<T>         m_costs_backup;
-    T m_converted_harris_eps;
+    T                 m_converted_harris_eps;
-    unsigned m_inf_row_index_for_tableau;
+    unsigned          m_inf_row_index_for_tableau;
-    bool m_bland_mode_tableau;
+    bool              m_bland_mode_tableau;
-    int_set m_left_basis_tableau;
+    int_set           m_left_basis_tableau;
-    unsigned m_bland_mode_threshold = 1000;
+    unsigned          m_bland_mode_threshold;
-    unsigned m_left_basis_repeated;
+    unsigned          m_left_basis_repeated;
-    vector<unsigned> m_leaving_candidates;
+    vector<unsigned>  m_leaving_candidates;
     //    T m_converted_harris_eps = convert_struct<T, double>::convert(this->m_settings.harris_feasibility_tolerance);
     std::list<unsigned> m_non_basis_list;
     void sort_non_basis();
@@ -76,10 +76,10 @@ public:
     //             choices.clear();
     //             choices.push_back(i);
     //             len = row_len;
-    //             if (my_random() % 10) break;
+    //             if (m_settings.random_next() % 10) break;
     //         } else if (row_len == len) {
     //             choices.push_back(i);
-    //             if (my_random() % 10) break;
+    //             if (m_settings.random_next() % 10) break;
     //         }
     //     }
 
@@ -89,7 +89,7 @@ public:
     //     if (choices.size() == 1)
     //         return choices[0];
         
-    //     unsigned k = my_random() % choices.size();
+    //     unsigned k = this->m_settings.random_next() % choices.size();
     //     return choices[k];
     //     #endif
     // }
@@ -287,7 +287,7 @@ public:
                 choices.clear();
                 choices.push_back(&rc);
             } else if (damage == num_of_non_free_basics &&
-                       this->m_A.m_columns[j].size() <= len && (my_random() % 2)) {
+                       this->m_A.m_columns[j].size() <= len && (this->m_settings.random_next() % 2)) {
                 choices.push_back(&rc);
                 len = this->m_A.m_columns[j].size();
             }
@@ -299,7 +299,7 @@ public:
             return -1;
         }
         const row_cell<T>* rc = choices.size() == 1? choices[0] :
-            choices[my_random() % choices.size()];
+            choices[this->m_settings.random_next() % choices.size()];
 
         a_ent = rc->m_value;
         return rc->m_j;
@@ -423,7 +423,7 @@ public:
 
     void find_feasible_solution();
 
-    bool is_tiny() const {return this->m_m < 10 && this->m_n < 20;}
+    // bool is_tiny() const {return this->m_m < 10 && this->m_n < 20;}
 
     void one_iteration();
     void one_iteration_tableau();
@@ -905,7 +905,9 @@ public:
                                   column_type_array,
                                   low_bound_values,
                                   upper_bound_values),
-        m_beta(A.row_count()) {
+        m_beta(A.row_count()),
+        m_epsilon_of_reduced_cost(T(1)/T(10000000)),
+        m_bland_mode_threshold(1000) {
 
         if (!(numeric_traits<T>::precise())) {
             m_converted_harris_eps = convert_struct<T, double>::convert(this->m_settings.harris_feasibility_tolerance);

src/util/lp/lp_primal_core_solver.hpp

@@ -199,7 +199,7 @@ int lp_primal_core_solver<T, X>::choose_entering_column_presize(unsigned number_
             entering_iter = non_basis_iter;
             if (number_of_benefitial_columns_to_go_over)
                 number_of_benefitial_columns_to_go_over--;
-        } else if (t == j_nz && my_random() % 2 == 0) {
+        } else if (t == j_nz && this->m_settings.random_next() % 2 == 0) {
             entering_iter = non_basis_iter;
         }
     }// while (number_of_benefitial_columns_to_go_over && initial_offset_in_non_basis != offset_in_nb);
@@ -268,7 +268,7 @@ template <typename T, typename X> int lp_primal_core_solver<T, X>::advance_on_so
         if (slope_at_entering * m_sign_of_entering_delta > - m_epsilon_of_reduced_cost) { // the slope started to increase infeasibility
             break;
         } else {
-            if ((numeric_traits<T>::precise() == false) || ( numeric_traits<T>::is_zero(slope_at_entering) && my_random() % 2 == 0)) {
+            if ((numeric_traits<T>::precise() == false) || ( numeric_traits<T>::is_zero(slope_at_entering) && this->m_settings.random_next() % 2 == 0)) {
                 // it is not cost benefitial to advance the delta more, so just break to increas the randomness
                 break;
             }
@@ -307,7 +307,7 @@ find_leaving_on_harris_theta(X const & harris_theta, X & t) {
     // we also know that harris_theta is limited, so we will find a leaving
     zero_harris_eps();
     unsigned steps = this->m_ed.m_index.size();
-    unsigned k = my_random() % steps;
+    unsigned k = this->m_settings.random_next() % steps;
     unsigned initial_k = k;
     do {
         unsigned i = this->m_ed.m_index[k];
@@ -398,7 +398,7 @@ template <typename T, typename X> int lp_primal_core_solver<T, X>::find_leaving_
         return find_leaving_and_t_with_breakpoints(entering, t);
     bool unlimited = true;
     unsigned steps = this->m_ed.m_index.size();
-    unsigned k = my_random() % steps;
+    unsigned k = this->m_settings.random_next() % steps;
     unsigned initial_k = k;
     unsigned row_min_nz = this->m_n() + 1;
     m_leaving_candidates.clear();
@@ -454,7 +454,7 @@ template <typename T, typename X> int lp_primal_core_solver<T, X>::find_leaving_
         t = ratio;
         return entering;
     }
-    k = my_random() % m_leaving_candidates.size();
+    k = this->m_settings.random_next() % m_leaving_candidates.size();
     return m_leaving_candidates[k];
 }
 
@@ -628,7 +628,7 @@ template <typename T, typename X>    void lp_primal_core_solver<T, X>::backup_an
 template <typename T, typename X>    void lp_primal_core_solver<T, X>::init_run() {
     this->m_basis_sort_counter = 0; // to initiate the sort of the basis
     this->set_total_iterations(0);
-    this->m_iters_with_no_cost_growing = 0;
+    this->iters_with_no_cost_growing() = 0;
     init_inf_set();
     if (this->current_x_is_feasible() && this->m_look_for_feasible_solution_only)
         return;
@@ -664,7 +664,7 @@ void lp_primal_core_solver<T, X>::advance_on_entering_equal_leaving(int entering
         this->init_lu();
         if (!this->find_x_by_solving()) {
             this->restore_x(entering, t * m_sign_of_entering_delta);
-            this->m_iters_with_no_cost_growing++;
+            this->iters_with_no_cost_growing()++;
             LP_OUT(this->m_settings, "failing in advance_on_entering_equal_leaving for entering = " << entering << std::endl);
             return;
         }
@@ -679,7 +679,7 @@ void lp_primal_core_solver<T, X>::advance_on_entering_equal_leaving(int entering
     if (need_to_switch_costs() ||!this->current_x_is_feasible()) {
         init_reduced_costs();
     }
-    this->m_iters_with_no_cost_growing = 0;
+    this->iters_with_no_cost_growing() = 0;
 }
 
 template <typename T, typename X>void lp_primal_core_solver<T, X>::advance_on_entering_and_leaving(int entering, int leaving, X & t) {
@@ -699,14 +699,14 @@ template <typename T, typename X>void lp_primal_core_solver<T, X>::advance_on_en
     if (!pivot_compare_result){;}
     else if (pivot_compare_result == 2) { // the sign is changed, cannot continue
         this->set_status(UNSTABLE);
-        this->m_iters_with_no_cost_growing++;
+        this->iters_with_no_cost_growing()++;
         return;
     } else {
         lean_assert(pivot_compare_result == 1);
         this->init_lu();
         if (this->m_factorization == nullptr || this->m_factorization->get_status() != LU_status::OK) {
             this->set_status(UNSTABLE);
-            this->m_iters_with_no_cost_growing++;
+            this->iters_with_no_cost_growing()++;
             return;
         }
     }
@@ -728,7 +728,7 @@ template <typename T, typename X>void lp_primal_core_solver<T, X>::advance_on_en
     }
 
     if (!is_zero(t)) {
-        this->m_iters_with_no_cost_growing = 0;
+        this->iters_with_no_cost_growing() = 0;
         init_infeasibility_after_update_x_if_inf(leaving);
     }
 
@@ -783,7 +783,7 @@ template <typename T, typename X> void lp_primal_core_solver<T, X>::advance_on_e
         this->init_lu();
         init_reduced_costs();
         if (refresh_result == 2) {
-            this->m_iters_with_no_cost_growing++;
+            this->iters_with_no_cost_growing()++;
             return;
         }
     }
@@ -833,7 +833,7 @@ template <typename T, typename X>  unsigned lp_primal_core_solver<T, X>::get_num
     if (ret == 0) {
         return 0;
     }
-    return std::max(static_cast<unsigned>(my_random() % ret), 1u);
+    return std::max(static_cast<unsigned>(this->m_settings.random_next() % ret), 1u);
 }
 
 template <typename T, typename X> void lp_primal_core_solver<T, X>::print_column_norms(std::ostream & out) {
@@ -934,7 +934,7 @@ template <typename T, typename X> unsigned lp_primal_core_solver<T, X>::solve()
              &&
              this->get_status() != INFEASIBLE
              &&
-             this->m_iters_with_no_cost_growing <= this->m_settings.max_number_of_iterations_with_no_improvements
+             this->iters_with_no_cost_growing() <= this->m_settings.max_number_of_iterations_with_no_improvements
              &&
              this->total_iterations() <= this->m_settings.max_total_number_of_iterations
              &&
@@ -961,7 +961,7 @@ template <typename T, typename X> void lp_primal_core_solver<T, X>::init_column_
     for (unsigned j = 0; j < this->m_n(); j++) {
         this->m_column_norms[j] = T(static_cast<int>(this->m_A.m_columns[j].size() + 1)) 
             
-            + T(static_cast<int>(my_random() % 10000)) / T(100000);
+            + T(static_cast<int>(this->m_settings.random_next() % 10000)) / T(100000);
     }
 }
 

src/util/lp/lp_primal_core_solver_instances.cpp

@@ -11,6 +11,7 @@
 #include "util/lp/lp_primal_core_solver.hpp"
 #include "util/lp/lp_primal_core_solver_tableau.hpp"
 namespace lean {
+
 template void lp_primal_core_solver<double, double>::find_feasible_solution();
 template void lean::lp_primal_core_solver<lean::mpq, lean::numeric_pair<lean::mpq> >::find_feasible_solution();
 
@@ -22,4 +23,5 @@ template void lean::lp_primal_core_solver<double, double>::clear_breakpoints();
 template bool lean::lp_primal_core_solver<lean::mpq, lean::mpq>::update_basis_and_x_tableau(int, int, lean::mpq const&);
 template bool lean::lp_primal_core_solver<double, double>::update_basis_and_x_tableau(int, int, double const&);
 template bool lean::lp_primal_core_solver<lean::mpq, lean::numeric_pair<lean::mpq> >::update_basis_and_x_tableau(int, int, lean::numeric_pair<lean::mpq> const&);
+
 }

src/util/lp/lp_primal_core_solver_tableau.hpp

@@ -62,7 +62,7 @@ template <typename T, typename X> int lp_primal_core_solver<T, X>::choose_enteri
             if (number_of_benefitial_columns_to_go_over)
                 number_of_benefitial_columns_to_go_over--;
         }
-        else if (t == j_nz && my_random() % 2 == 0) {
+        else if (t == j_nz && this->m_settings.random_next() % 2 == 0) {
             entering_iter = non_basis_iter;
         }
     }// while (number_of_benefitial_columns_to_go_over && initial_offset_in_non_basis != offset_in_nb);
@@ -169,7 +169,7 @@ unsigned lp_primal_core_solver<T, X>::solve_with_tableau() {
              &&
              this->get_status() != INFEASIBLE
              &&
-             this->m_iters_with_no_cost_growing <= this->m_settings.max_number_of_iterations_with_no_improvements
+             this->iters_with_no_cost_growing() <= this->m_settings.max_number_of_iterations_with_no_improvements
              &&
              this->total_iterations() <= this->m_settings.max_total_number_of_iterations
              &&
@@ -202,7 +202,7 @@ template <typename T, typename X>void lp_primal_core_solver<T, X>::advance_on_en
         }
         this->update_basis_and_x_tableau(entering, leaving, t);
         lean_assert(this->A_mult_x_is_off() == false);
-        this->m_iters_with_no_cost_growing = 0;
+        this->iters_with_no_cost_growing() = 0;
     } else {
         this->pivot_column_tableau(entering, this->m_basis_heading[leaving]);
         this->change_basis(entering, leaving);
@@ -233,7 +233,7 @@ void lp_primal_core_solver<T, X>::advance_on_entering_equal_leaving_tableau(int
     if (need_to_switch_costs()) {
         init_reduced_costs_tableau();
     }
-    this->m_iters_with_no_cost_growing = 0;
+    this->iters_with_no_cost_growing() = 0;
 }
 template <typename T, typename X> int lp_primal_core_solver<T, X>::find_leaving_and_t_tableau(unsigned entering, X & t) {
     unsigned k = 0;
@@ -293,7 +293,7 @@ template <typename T, typename X> int lp_primal_core_solver<T, X>::find_leaving_
     }
     if (m_leaving_candidates.size() == 1)
         return m_leaving_candidates[0];
-    k = my_random() % m_leaving_candidates.size();
+    k = this->m_settings.random_next() % m_leaving_candidates.size();
     return m_leaving_candidates[k];
 }
 template <typename T, typename X> void lp_primal_core_solver<T, X>::init_run_tableau() {
@@ -302,7 +302,7 @@ template <typename T, typename X> void lp_primal_core_solver<T, X>::init_run_tab
         lean_assert(basis_columns_are_set_correctly());
         this->m_basis_sort_counter = 0; // to initiate the sort of the basis
         this->set_total_iterations(0);
-        this->m_iters_with_no_cost_growing = 0;
+        this->iters_with_no_cost_growing() = 0;
 		lean_assert(this->inf_set_is_correct());
         if (this->current_x_is_feasible() && this->m_look_for_feasible_solution_only)
             return;
@@ -315,7 +315,7 @@ template <typename T, typename X> void lp_primal_core_solver<T, X>::init_run_tab
             this->m_column_norm_update_counter = 0;
             init_column_norms();
         }
-        if (this->m_settings.m_simplex_strategy == simplex_strategy_enum::tableau_rows)
+        if (this->m_settings.simplex_strategy() == simplex_strategy_enum::tableau_rows)
             init_tableau_rows();
         lean_assert(this->reduced_costs_are_correct_tableau());
         lean_assert(!this->need_to_pivot_to_basis_tableau());

src/util/lp/lp_primal_simplex.h

@@ -15,7 +15,7 @@
 namespace lean {
 template <typename T, typename X>
 class lp_primal_simplex: public lp_solver<T, X> {
-    lp_primal_core_solver<T, X> * m_core_solver = nullptr;
+    lp_primal_core_solver<T, X> * m_core_solver;
     vector<X> m_low_bounds;
 private:
     unsigned original_rows() { return this->m_external_rows_to_core_solver_rows.size(); }
@@ -28,7 +28,7 @@ private:
 
     void set_scaled_costs();
 public:
-    lp_primal_simplex() {}
+    lp_primal_simplex(): m_core_solver(nullptr) {}
 
     column_info<T> * get_or_create_column_info(unsigned column);
 

src/util/lp/lp_settings.h

@@ -25,9 +25,10 @@ enum class column_type  {
         };
 
 enum class simplex_strategy_enum {
+    undecided = 3,
     tableau_rows = 0,
     tableau_costs = 1,
-    no_tableau = 2
+    lu = 2
 };
 
 std::string column_type_to_string(column_type t);
@@ -70,8 +71,6 @@ template <typename X> bool is_epsilon_small(const X & v, const double& eps);
 
 int get_millisecond_count();
 int get_millisecond_span(int start_time);
-unsigned my_random();
-void my_random_init(long unsigned seed);
 
 
 class lp_resource_limit {
@@ -105,49 +104,50 @@ private:
     default_lp_resource_limit m_default_resource_limit;
     lp_resource_limit* m_resource_limit;
     // used for debug output
-    std::ostream* m_debug_out = &std::cout;
+    std::ostream* m_debug_out;
     // used for messages, for example, the computation progress messages
-    std::ostream* m_message_out = &std::cout;
+    std::ostream* m_message_out;
 
     stats  m_stats;
+    random_gen m_rand;
 
 public:
-    unsigned reps_in_scaler = 20;
+    unsigned reps_in_scaler;
     // when the absolute value of an element is less than pivot_epsilon
     // in pivoting, we treat it as a zero
-    double pivot_epsilon = 0.00000001;
+    double pivot_epsilon;
     // see Chatal, page 115
-    double positive_price_epsilon = 1e-7;
+    double positive_price_epsilon;
     // a quatation "if some choice of the entering vairable leads to an eta matrix
     // whose diagonal element in the eta column is less than e2 (entering_diag_epsilon) in magnitude, the this choice is rejected ...
-    double entering_diag_epsilon = 1e-8;
+    double entering_diag_epsilon;
-    int c_partial_pivoting = 10; // this is the constant c from page 410
+    int c_partial_pivoting; // this is the constant c from page 410
-    unsigned depth_of_rook_search = 4;
+    unsigned depth_of_rook_search;
-    bool using_partial_pivoting = true;
+    bool using_partial_pivoting;
     // dissertation of Achim Koberstein
     // if Bx - b is different at any component more that refactor_epsilon then we refactor
-    double refactor_tolerance = 1e-4;
+    double refactor_tolerance;
-    double pivot_tolerance = 1e-6;
+    double pivot_tolerance;
-    double zero_tolerance = 1e-12;
+    double zero_tolerance;
-    double drop_tolerance = 1e-14;
+    double drop_tolerance;
-    double tolerance_for_artificials = 1e-4;
+    double tolerance_for_artificials;
-    double can_be_taken_to_basis_tolerance = 0.00001;
+    double can_be_taken_to_basis_tolerance;
 
-    unsigned percent_of_entering_to_check = 5; // we try to find a profitable column in a percentage of the columns
+    unsigned percent_of_entering_to_check; // we try to find a profitable column in a percentage of the columns
-    bool use_scaling = true;
+    bool use_scaling;
-    double scaling_maximum = 1;
+    double scaling_maximum;
-    double scaling_minimum = 0.5;
+    double scaling_minimum;
-    double harris_feasibility_tolerance = 1e-7; // page 179 of Istvan Maros
+    double harris_feasibility_tolerance; // page 179 of Istvan Maros
-    double ignore_epsilon_of_harris = 10e-5;
+    double ignore_epsilon_of_harris;
-    unsigned max_number_of_iterations_with_no_improvements = 2000000;
+    unsigned max_number_of_iterations_with_no_improvements;
-    unsigned max_total_number_of_iterations = 20000000;
+    unsigned max_total_number_of_iterations;
-    double time_limit = std::numeric_limits<double>::max(); // the maximum time limit of the total run time in seconds
+    double time_limit; // the maximum time limit of the total run time in seconds
     // dual section
-    double dual_feasibility_tolerance = 1e-7; // // page 71 of the PhD thesis of Achim Koberstein
+    double dual_feasibility_tolerance; // // page 71 of the PhD thesis of Achim Koberstein
-    double primal_feasibility_tolerance = 1e-7; // page 71 of the PhD thesis of Achim Koberstein
+    double primal_feasibility_tolerance; // page 71 of the PhD thesis of Achim Koberstein
-    double relative_primal_feasibility_tolerance = 1e-9; // page 71 of the PhD thesis of Achim Koberstein
+    double relative_primal_feasibility_tolerance; // page 71 of the PhD thesis of Achim Koberstein
 
-    bool m_bound_propagation = true;
+    bool m_bound_propagation;
     
     bool bound_progation() const {
         return m_bound_propagation;
@@ -157,7 +157,52 @@ public:
         return m_bound_propagation;
     }
     
-    lp_settings() : m_default_resource_limit(*this), m_resource_limit(&m_default_resource_limit) {}
+    lp_settings() : m_default_resource_limit(*this),
+                    m_resource_limit(&m_default_resource_limit),
+                    m_debug_out( &std::cout),
+                    m_message_out(&std::cout),
+                    reps_in_scaler(20),
+                    pivot_epsilon(0.00000001),
+                    positive_price_epsilon(1e-7),
+                    entering_diag_epsilon ( 1e-8),
+                    c_partial_pivoting ( 10), // this is the constant c from page 410
+                    depth_of_rook_search ( 4),
+                    using_partial_pivoting ( true),
+                    // dissertation of Achim Koberstein
+                    // if Bx - b is different at any component more that refactor_epsilon then we refactor
+                    refactor_tolerance ( 1e-4),
+                    pivot_tolerance ( 1e-6),
+                    zero_tolerance ( 1e-12),
+                    drop_tolerance ( 1e-14),
+                    tolerance_for_artificials ( 1e-4),
+                    can_be_taken_to_basis_tolerance ( 0.00001),
+                  
+                    percent_of_entering_to_check ( 5),// we try to find a profitable column in a percentage of the columns
+                    use_scaling ( true),
+                    scaling_maximum ( 1),
+                    scaling_minimum ( 0.5),
+                    harris_feasibility_tolerance ( 1e-7), // page 179 of Istvan Maros
+                    ignore_epsilon_of_harris ( 10e-5),
+                    max_number_of_iterations_with_no_improvements ( 2000000),
+                    max_total_number_of_iterations ( 20000000),
+                    time_limit ( std::numeric_limits<double>::max()), // the maximum time limit of the total run time in seconds
+                    // dual section
+                    dual_feasibility_tolerance ( 1e-7), // // page 71 of the PhD thesis of Achim Koberstein
+                    primal_feasibility_tolerance ( 1e-7), // page 71 of the PhD thesis of Achim Koberstein
+                    relative_primal_feasibility_tolerance ( 1e-9), // page 71 of the PhD thesis of Achim Koberstein
+                    m_bound_propagation ( true),
+                    presolve_with_double_solver_for_lar(true),
+                    m_simplex_strategy(simplex_strategy_enum::tableau_rows),
+                    report_frequency(1000),
+                    print_statistics(false),
+                    column_norms_update_frequency(12000),
+                    scale_with_ratio(true),
+                    density_threshold(0.7),
+                    use_breakpoints_in_feasibility_search(false),
+                    max_row_length_for_bound_propagation(300),
+                    backup_costs(true),
+                    column_number_threshold_for_using_lu_in_lar_solver(4000)
+    {}
 
     void set_resource_limit(lp_resource_limit& lim) { m_resource_limit = &lim; }
     bool get_cancel_flag() const { return m_resource_limit->get_cancel_flag(); }
@@ -226,8 +271,8 @@ public:
         return is_eps_small_general<T>(t, tolerance_for_artificials);
     }
     // the method of lar solver to use
-    bool presolve_with_double_solver_for_lar = true;
+    bool presolve_with_double_solver_for_lar;
-    simplex_strategy_enum m_simplex_strategy = simplex_strategy_enum::tableau_rows;
+    simplex_strategy_enum m_simplex_strategy;
     simplex_strategy_enum simplex_strategy() const {
         return m_simplex_strategy;
     }
@@ -236,27 +281,33 @@ public:
         return m_simplex_strategy;
     }
 
+	bool use_lu() const {
+		return m_simplex_strategy == simplex_strategy_enum::lu;
+	}
+
     bool use_tableau() const {
-        return m_simplex_strategy != simplex_strategy_enum::no_tableau;
+		return m_simplex_strategy == simplex_strategy_enum::tableau_rows ||
+			m_simplex_strategy == simplex_strategy_enum::tableau_costs;
     }
 
     bool use_tableau_rows() const {
         return m_simplex_strategy == simplex_strategy_enum::tableau_rows;
     }
     
-    int report_frequency = 1000;
+    int report_frequency;
-    bool print_statistics = false;
+    bool print_statistics;
-    unsigned column_norms_update_frequency = 12000;
+    unsigned column_norms_update_frequency;
-    bool scale_with_ratio = true;
+    bool scale_with_ratio;
-    double density_threshold = 0.7; // need to tune it up, todo
+    double density_threshold; // need to tune it up, todo
 #ifdef LEAN_DEBUG
     static unsigned ddd; // used for debugging    
 #endif
-    bool use_breakpoints_in_feasibility_search = false;
+    bool use_breakpoints_in_feasibility_search;
-    unsigned random_seed = 1;
+    unsigned random_next() { return m_rand(); }
-    static unsigned long random_next;
+    void random_seed(unsigned s) { m_rand.set_seed(s); }
-    unsigned max_row_length_for_bound_propagation = 300;
+    unsigned max_row_length_for_bound_propagation;
-    bool backup_costs = true;
+    bool backup_costs;
+    unsigned column_number_threshold_for_using_lu_in_lar_solver;
 }; // end of lp_settings class
 
 

src/util/lp/lp_settings.hpp

@@ -67,15 +67,6 @@ int get_millisecond_span(int start_time) {
 
 
 
-void my_random_init(long unsigned seed) {
-    lp_settings::random_next = seed;   
-}
-
-unsigned my_random() {
-    lp_settings::random_next = lp_settings::random_next * 1103515245 + 12345;
-    return((unsigned)(lp_settings::random_next/65536) % 32768);
-}
-
 template <typename T>
 bool vectors_are_equal(T * a, vector<T>  &b, unsigned n) {
     if (numeric_traits<T>::precise()) {
@@ -126,7 +117,6 @@ bool vectors_are_equal(const vector<T> & a, const vector<T>  &b) {
     }
     return true;
 }
-unsigned long lp_settings::random_next = 1;
 #ifdef LEAN_DEBUG
 unsigned lp_settings::ddd = 0;
 #endif

src/util/lp/lp_solver.h

@@ -39,10 +39,10 @@ protected:
     T get_column_cost_value(unsigned j, column_info<T> * ci) const;
 public:
     unsigned m_total_iterations;
-    static_matrix<T, X>* m_A = nullptr; // this is the matrix of constraints
+    static_matrix<T, X>* m_A; // this is the matrix of constraints
     vector<T> m_b; // the right side vector
-    unsigned m_first_stage_iterations = 0;
+    unsigned m_first_stage_iterations;
-    unsigned m_second_stage_iterations = 0;
+    unsigned m_second_stage_iterations;
     std::unordered_map<unsigned, lp_constraint<T, X>> m_constraints;
     std::unordered_map<var_index, column_info<T>*> m_map_from_var_index_to_column_info;
     std::unordered_map<unsigned, std::unordered_map<unsigned, T> > m_A_values;
@@ -52,8 +52,8 @@ public:
     std::unordered_map<unsigned, unsigned> m_core_solver_columns_to_external_columns;
     vector<T> m_column_scale;
     std::unordered_map<unsigned, std::string>  m_name_map;
-    unsigned m_artificials = 0;
+    unsigned m_artificials;
-    unsigned m_slacks = 0;
+    unsigned m_slacks;
     vector<column_type> m_column_types;
     vector<T> m_costs;
     vector<T> m_x;
@@ -63,10 +63,17 @@ public:
     vector<int> m_heading;
 
 
-    lp_status m_status = lp_status::UNKNOWN;
+    lp_status m_status;
 
     lp_settings m_settings;
-    lp_solver() {}
+    lp_solver():
+        m_A(nullptr), // this is the matrix of constraints
+        m_first_stage_iterations (0),
+        m_second_stage_iterations (0),
+        m_artificials (0),
+        m_slacks (0),
+        m_status(lp_status::UNKNOWN)
+    {}
 
     unsigned row_count() const { return this->m_A->row_count(); }
 
@@ -232,14 +239,6 @@ protected:
         out << "extended A[" << this->m_A->row_count() << "," << this->m_A->column_count() << "]" << std::endl;
     }
 
-    struct row_tighten_stats {
-        unsigned n_of_new_bounds = 0;
-        unsigned n_of_fixed = 0;
-        bool is_obsolete = false;
-    };
-    
-
-
 public:
     lp_settings & settings() { return m_settings;}
     void print_model(std::ostream & s) const {

src/util/lp/lu.h

@@ -110,25 +110,25 @@ enum class LU_status { OK, Degenerated};
 // Using Suhl-Suhl method described in the dissertation of Achim Koberstein, Chapter 5
 template <typename T, typename X>
 class lu {
-    LU_status m_status = LU_status::OK;
+    LU_status m_status;
 public:
     // the fields
-    unsigned m_dim;
+    unsigned                   m_dim;
     static_matrix<T, X> const &m_A;
-    permutation_matrix<T, X> m_Q;
+    permutation_matrix<T, X>   m_Q;
-    permutation_matrix<T, X> m_R;
+    permutation_matrix<T, X>   m_R;
-    permutation_matrix<T, X> m_r_wave;
+    permutation_matrix<T, X>   m_r_wave;
-    sparse_matrix<T, X> m_U;
+    sparse_matrix<T, X>        m_U;
     square_dense_submatrix<T, X>* m_dense_LU;
     
     vector<tail_matrix<T, X> *> m_tail;
-    lp_settings & m_settings;
+    lp_settings &               m_settings;
-    bool m_failure = false;
+    bool                        m_failure;
-    indexed_vector<T> m_row_eta_work_vector;
+    indexed_vector<T>           m_row_eta_work_vector;
-    indexed_vector<T> m_w_for_extension;
+    indexed_vector<T>           m_w_for_extension;
-    indexed_vector<T> m_y_copy;
+    indexed_vector<T>           m_y_copy;
-    indexed_vector<unsigned> m_ii; //to optimize the work with the m_index fields
+    indexed_vector<unsigned>    m_ii; //to optimize the work with the m_index fields
-    unsigned m_refactor_counter = 0;
+    unsigned                    m_refactor_counter;
     // constructor
     // if A is an m by n matrix then basis has length m and values in [0,n); the values are all different
     // they represent the set of m columns

src/util/lp/lu.hpp

@@ -111,6 +111,7 @@ template <typename T, typename X>
 lu<T, X>::lu(static_matrix<T, X> const & A,
              vector<unsigned>& basis,
              lp_settings & settings):
+    m_status(LU_status::OK),
     m_dim(A.row_count()),
     m_A(A),
     m_Q(m_dim),
@@ -118,7 +119,9 @@ lu<T, X>::lu(static_matrix<T, X> const & A,
     m_r_wave(m_dim),
     m_U(A, basis), // create the square matrix that eventually will be factorized
     m_settings(settings),
-    m_row_eta_work_vector(A.row_count()){
+    m_failure(false),
+    m_row_eta_work_vector(A.row_count()),
+    m_refactor_counter(0) {
     lean_assert(!(numeric_traits<T>::precise() && settings.use_tableau()));
 #ifdef LEAN_DEBUG
     debug_test_of_basis(A, basis);
@@ -602,13 +605,13 @@ void lu<T, X>::process_column(int j) {
     unsigned pi, pj;
     bool success = m_U.get_pivot_for_column(pi, pj, m_settings.c_partial_pivoting, j);
     if (!success) {
-        LP_OUT(m_settings, "get_pivot returned false: cannot find the pivot for column " << j << std::endl);
+        //        LP_OUT(m_settings, "get_pivot returned false: cannot find the pivot for column " << j << std::endl);
         m_failure = true;
         return;
     }
 
     if (static_cast<int>(pi) == -1) {
-        LP_OUT(m_settings, "cannot find the pivot for column " << j << std::endl);
+        // LP_OUT(m_settings, "cannot find the pivot for column " << j << std::endl);
         m_failure = true;
         return;
     }

src/util/lp/mps_reader.h

@@ -93,22 +93,28 @@ template <typename T, typename X>
 class mps_reader {
     enum row_type { Cost, Less_or_equal, Greater_or_equal, Equal };
     struct bound {
-        bool m_low_is_set = true;
+        T    m_low;
-        T m_low;
+        T    m_upper;
-        bool m_upper_is_set = false;
+        bool m_low_is_set;
-        T m_upper;
+        bool m_upper_is_set;
-        bool m_value_is_fixed = false;
+        bool m_value_is_fixed;
-        T m_fixed_value;
+        T    m_fixed_value;
-        bool m_free = false;
+        bool m_free;
         // constructor
-        bound() : m_low(numeric_traits<T>::zero()) {} // it seems all mps files I have seen have the default low value 0 on a variable
+        bound() : m_low(numeric_traits<T>::zero()),
+                  m_low_is_set(true),
+                  m_upper_is_set(false),
+                  m_value_is_fixed(false),
+                  m_free(false) {} // it seems all mps files I have seen have the default low value 0 on a variable
     };
 
     struct column {
         std::string m_name;
-        bound * m_bound = nullptr;
+        bound * m_bound;
         unsigned m_index;
-        column(std::string name, unsigned index): m_name(name), m_index(index) {
+        column(std::string name, unsigned index): m_name(name),
+                                                  m_bound(nullptr),
+                                                  m_index(index) {
         }
     };
 
@@ -116,15 +122,21 @@ class mps_reader {
         row_type m_type;
         std::string m_name;
         std::unordered_map<std::string, T> m_row_columns;
-        T m_right_side = numeric_traits<T>::zero();
         unsigned m_index;
-        T m_range = numeric_traits<T>::zero();
+        T m_right_side;
-        row(row_type type, std::string name, unsigned index) : m_type(type), m_name(name), m_index(index) {
+        T m_range;
+        row(row_type type, std::string name, unsigned index) :
+            m_type(type),
+            m_name(name),
+            m_index(index),
+            m_right_side(zero_of_type<T>()),
+            m_range(zero_of_type<T>())
+        {
         }
     };
 
+    bool m_is_OK;
     std::string m_file_name;
-    bool m_is_OK = true;
     std::unordered_map<std::string, row *> m_rows;
     std::unordered_map<std::string, column *> m_columns;
     std::unordered_map<std::string, unsigned> m_names_to_var_index;
@@ -133,9 +145,9 @@ class mps_reader {
     std::string m_cost_row_name;
     std::ifstream m_file_stream;
     // needed to adjust the index row
-    unsigned m_cost_line_count = 0;
+    unsigned m_cost_line_count;
-    unsigned m_line_number = 0;
+    unsigned m_line_number;
-    std::ostream * m_message_stream = & std::cout;
+    std::ostream * m_message_stream;
 
     void set_m_ok_to_false() {
         *m_message_stream << "setting m_is_OK to false" << std::endl;
@@ -737,8 +749,12 @@ public:
     }
 
     mps_reader(std::string file_name):
-        m_file_name(file_name), m_file_stream(file_name) {
+        m_is_OK(true),
-    }
+        m_file_name(file_name), 
+        m_file_stream(file_name),
+        m_cost_line_count(0),
+        m_line_number(0),
+        m_message_stream(& std::cout) {}
     void read() {
         if (!m_file_stream.is_open()){
             set_m_ok_to_false();
@@ -784,7 +800,7 @@ public:
         auto it = m_names_to_var_index.find(s);
         if (it != m_names_to_var_index.end())
             return it->second;
-        unsigned ret = m_names_to_var_index.size();
+        unsigned ret = static_cast<unsigned>(m_names_to_var_index.size());
         m_names_to_var_index[s] = ret;
         return ret;
     }

src/util/lp/numeric_pair.h

@@ -101,16 +101,14 @@ struct numeric_pair {
 
     numeric_pair(T xp, T yp) : x(xp), y(yp) {}
 
-
     template <typename X>
     numeric_pair(const X & n) : x(n), y(0) {
     }
     
-    template <typename X>
+    numeric_pair(const numeric_pair<T> & n) : x(n.x), y(n.y) {}
-    numeric_pair(const numeric_pair<X> & n) : x(n.x), y(n.y) {}
     
     template <typename X, typename Y>
-    numeric_pair(X xp, Y yp) : numeric_pair(convert_struct<T, X>::convert(xp), convert_struct<T, Y>::convert(yp)) {}
+    numeric_pair(X xp, Y yp) : x(convert_struct<T, X>::convert(xp)), y(convert_struct<T, Y>::convert(yp)) {}
 
     bool operator<(const numeric_pair& a) const {
         return x < a.x || (x == a.x && y < a.y);

src/util/lp/permutation_matrix.h

@@ -132,42 +132,4 @@ class permutation_matrix : public tail_matrix<T, X> {
         
     }; // end of the permutation class
 
-#ifdef LEAN_DEBUG
-template <typename T, typename X>
-class permutation_generator {
-    unsigned m_n;
-    permutation_generator* m_lower;
-    bool m_done = false;
-    permutation_matrix<T, X> m_current;
-    unsigned m_last;
-public:
-    permutation_generator(unsigned n);
-    permutation_generator(const permutation_generator & o);
-    bool move_next();
-
-    ~permutation_generator() {
-        if (m_lower != nullptr) {
-            delete m_lower;
-        }
-    }
-
-    permutation_matrix<T, X> *current() {
-        return &m_current;
-    }
-};
-
-template <typename T, typename X>
-inline unsigned number_of_inversions(permutation_matrix<T, X> & p);
-
-template <typename T, typename X>
-int sign(permutation_matrix<T, X> & p) {
-    return is_even(number_of_inversions(p))? 1: -1;
-}
-
-template <typename T, typename X>
-T det_val_on_perm(permutation_matrix<T, X>* u, const matrix<T, X>& m);
-
-template <typename T, typename X>
-T determinant(const matrix<T, X>& m);
-#endif
 }

src/util/lp/permutation_matrix.hpp

@@ -320,100 +320,4 @@ template <typename T, typename X> bool permutation_matrix<T, X>::is_identity() c
 }
 
 
-#ifdef LEAN_DEBUG
-template <typename T, typename X>
-permutation_generator<T, X>::permutation_generator(unsigned n): m_n(n), m_current(n) {
-    lean_assert(n > 0);
-    if (n > 1) {
-        m_lower = new permutation_generator(n - 1);
-    } else {
-        m_lower = nullptr;
-    }
-
-    m_last = 0;
-}
-
-template <typename T, typename X>
-permutation_generator<T, X>::permutation_generator(const permutation_generator & o): m_n(o.m_n), m_done(o.m_done), m_current(o.m_current), m_last(o.m_last) {
-    if (m_lower != nullptr) {
-        m_lower = new permutation_generator(o.m_lower);
-    } else {
-        m_lower = nullptr;
-    }
-}
-
-template <typename T, typename X> bool
-permutation_generator<T, X>::move_next() {
-    if (m_done) {
-        return false;
-    }
-
-    if (m_lower == nullptr) {
-        if (m_last == 0) {
-            m_last++;
-            return true;
-        } else {
-            m_done = true;
-            return false;
-        }
-    } else {
-        if (m_last < m_n && m_last > 0) {
-            m_current[m_last - 1] = m_current[m_last];
-            m_current[m_last] = m_n - 1;
-            m_last++;
-            return true;
-        } else {
-            if (m_lower -> move_next()) {
-                auto lower_curr = m_lower -> current();
-                for ( unsigned i = 1; i < m_n; i++ ){
-                    m_current[i] = (*lower_curr)[i - 1];
-                }
-                m_current[0] = m_n - 1;
-                m_last = 1;
-                return true;
-            } else {
-                m_done = true;
-                return false;
-            }
-        }
-    }
-}
-
-template <typename T, typename X>
-inline unsigned number_of_inversions(permutation_matrix<T, X> & p) {
-    unsigned ret = 0;
-    unsigned n = p.size();
-    for (unsigned i = 0; i < n; i++) {
-        for (unsigned j = i + 1; j < n; j++) {
-            if (p[i] > p[j]) {
-                ret++;
-            }
-        }
-    }
-    return ret;
-}
-
-template <typename T, typename X>
-T det_val_on_perm(permutation_matrix<T, X>* u, const matrix<T, X>& m) {
-    unsigned n = m.row_count();
-    T ret = numeric_traits<T>::one();
-    for (unsigned i = 0; i < n; i++) {
-        unsigned j = (*u)[i];
-        ret *=  m(i, j);
-    }
-    return ret * sign(*u);
-}
-
-template <typename T, typename X>
-T determinant(const matrix<T, X>& m) {
-    lean_assert(m.column_count() == m.row_count());
-    unsigned n = m.row_count();
-    permutation_generator<T, X> allp(n);
-    T ret = numeric_traits<T>::zero();
-    while (allp.move_next()){
-        ret += det_val_on_perm(allp.current(), m);
-    }
-    return ret;
-}
-#endif
 }

src/util/lp/permutation_matrix_instances.cpp

@@ -46,11 +46,6 @@ template void lean::permutation_matrix<lean::mpq, lean::numeric_pair<lean::mpq>
 template void lean::permutation_matrix<lean::mpq, lean::numeric_pair<lean::mpq> >::apply_reverse_from_left_to_T(vector<lean::mpq>&);
 template void lean::permutation_matrix<lean::mpq, lean::numeric_pair<lean::mpq> >::apply_reverse_from_right_to_T(vector<lean::mpq >&);
 template void lean::permutation_matrix<double, double>::multiply_by_permutation_from_right(lean::permutation_matrix<double, double>&);
-
-#ifdef LEAN_DEBUG
-template bool lean::permutation_generator<double, double>::move_next();
-template lean::permutation_generator<double, double>::permutation_generator(unsigned int);
-#endif
 template lean::permutation_matrix<double, double>::permutation_matrix(unsigned int);
 template void lean::permutation_matrix<double, double>::apply_reverse_from_left_to_X(vector<double> &);
 template void  lean::permutation_matrix< lean::mpq, lean::mpq>::apply_reverse_from_left_to_X(vector<lean::mpq> &);

src/util/lp/random_updater.h

@@ -16,12 +16,14 @@ namespace lean {
 template <typename T> struct numeric_pair; // forward definition
 class lar_core_solver; // forward definition
 class random_updater {
-    unsigned range = 100000;
     struct interval {
-        bool upper_bound_is_set = false;
+        bool upper_bound_is_set;
         numeric_pair<mpq> upper_bound;
-        bool low_bound_is_set = false;
+        bool low_bound_is_set;
         numeric_pair<mpq> low_bound;
+        interval() : upper_bound_is_set(false),
+                     low_bound_is_set(false) {}
+        
         void set_low_bound(const numeric_pair<mpq> & v) {
             if (low_bound_is_set) {
                 low_bound = std::max(v, low_bound);
@@ -58,6 +60,7 @@ class random_updater {
     };
     std::set<var_index> m_var_set;
     lar_core_solver & m_core_solver;
+    unsigned range;
     linear_combination_iterator<mpq>* m_column_j; // the actual column
     interval find_shift_interval(unsigned j);
     interval get_interval_of_non_basic_var(unsigned j);

src/util/lp/random_updater.hpp

@@ -12,7 +12,9 @@ namespace lean {
 
 random_updater::random_updater(
                                lar_core_solver & lar_core_solver,
-                               const vector<unsigned> & column_indices) : m_core_solver(lar_core_solver) {
+                               const vector<unsigned> & column_indices) :
+    m_core_solver(lar_core_solver),
+    range(100000) {
     for (unsigned j : column_indices)
         add_column_to_sets(j);
 }
@@ -134,7 +136,7 @@ void random_updater::shift_var(unsigned j, interval & r) {
 }
 
 numeric_pair<mpq> random_updater::get_random_from_interval(interval & r) {
-    unsigned rand = my_random();
+    unsigned rand = m_core_solver.settings().random_next();
     if ((!r.low_bound_is_set)  && (!r.upper_bound_is_set))
         return numeric_pair<mpq>(rand % range, 0);
     if (r.low_bound_is_set  && (!r.upper_bound_is_set))

src/util/lp/sparse_matrix.h

@@ -30,10 +30,10 @@ class sparse_matrix
 #endif
 {
     struct col_header {
-        unsigned m_shortened_markovitz = 0;
+        unsigned m_shortened_markovitz;
         vector<indexed_value<T>> m_values; // the actual column values
 
-        col_header()  {}
+        col_header(): m_shortened_markovitz(0) {}
 
         void shorten_markovich_by_one() {
             m_shortened_markovitz++;
@@ -44,17 +44,17 @@ class sparse_matrix
         }
     };
 
-    unsigned m_n_of_active_elems = 0;
+    unsigned                          m_n_of_active_elems;
     binary_heap_upair_queue<unsigned> m_pivot_queue;
 public:
     vector<vector<indexed_value<T>>>  m_rows;
-    vector<col_header> m_columns;
+    vector<col_header>                m_columns;
-    permutation_matrix<T, X>  m_row_permutation;
+    permutation_matrix<T, X>          m_row_permutation;
-    permutation_matrix<T, X>  m_column_permutation;
+    permutation_matrix<T, X>          m_column_permutation;
     // m_work_pivot_vector[j] = offset of elementh of j-th column in the row we are pivoting to
     // if the column is not present then m_work_pivot_vector[j] is -1
-    vector<int> m_work_pivot_vector;
+    vector<int>                       m_work_pivot_vector;
-    vector<bool> m_processed;
+    vector<bool>                      m_processed;
     unsigned get_n_of_active_elems() const { return m_n_of_active_elems; }
 
 #ifdef LEAN_DEBUG

src/util/lp/sparse_matrix.hpp

@@ -36,6 +36,7 @@ void sparse_matrix<T, X>::copy_B(static_matrix<T, X> const &A, vector<unsigned>
 // constructor that copies columns of the basis from A
 template <typename T, typename X>
 sparse_matrix<T, X>::sparse_matrix(static_matrix<T, X> const &A, vector<unsigned> & basis) :
+    m_n_of_active_elems(0),
     m_pivot_queue(A.row_count()),
     m_row_permutation(A.row_count()),
     m_column_permutation(A.row_count()),

src/util/lp/sparse_vector.h

@@ -20,7 +20,7 @@ public:
     }
 #ifdef LEAN_DEBUG
     T operator[] (unsigned i) const {
-        for (auto t : m_data) {
+        for (auto &t : m_data) {
             if (t.first == i) return t.second;
         }
         return numeric_traits<T>::zero();

src/util/lp/square_dense_submatrix.h

@@ -42,7 +42,7 @@ public:
     unsigned m_index_start;
     unsigned m_dim;
     vector<T> m_v;
-    sparse_matrix<T, X> * m_parent = nullptr;
+    sparse_matrix<T, X> * m_parent;
     permutation_matrix<T, X>  m_row_permutation;
     indexed_vector<T> m_work_vector;
 public:

src/util/lp/static_matrix.hpp

@@ -29,7 +29,6 @@ template <typename T, typename X> void static_matrix<T, X>::scan_row_ii_to_offse
 
 
 template <typename T, typename X> bool static_matrix<T, X>::pivot_row_to_row_given_cell(unsigned i, column_cell & c, unsigned pivot_col) {
-    //    std::cout << "ddd = " << ++lp_settings::ddd<< std::endl;
     unsigned ii = c.m_i;
     lean_assert(i < row_count() && ii < column_count());
     lean_assert(i != ii);

src/util/lp/ul_pair.h

@@ -32,14 +32,14 @@ inline bool compare(const std::pair<mpq, var_index> & a, const std::pair<mpq, va
 }
 
 class ul_pair {
-    constraint_index m_low_bound_witness = static_cast<constraint_index>(-1);
+    constraint_index m_low_bound_witness;
-    constraint_index m_upper_bound_witness = static_cast<constraint_index>(-1);
+    constraint_index m_upper_bound_witness;
 public:
     constraint_index& low_bound_witness() {return m_low_bound_witness;}
     constraint_index low_bound_witness() const {return m_low_bound_witness;}
     constraint_index& upper_bound_witness() { return m_upper_bound_witness;}
     constraint_index upper_bound_witness() const {return m_upper_bound_witness;}
-    row_index m_i = static_cast<row_index>(-1);
+    row_index m_i;
     bool operator!=(const ul_pair & p) const {
         return !(*this == p);
     }
@@ -50,8 +50,15 @@ public:
             m_i == p.m_i;
     }
 	// empty constructor
-	ul_pair(){}
+	ul_pair() :
-    ul_pair(row_index ri) : m_i(ri) {}
+        m_low_bound_witness(static_cast<constraint_index>(-1)),
+        m_upper_bound_witness(static_cast<constraint_index>(-1)),
+        m_i(static_cast<row_index>(-1))
+{}
+    ul_pair(row_index ri) :
+        m_low_bound_witness(static_cast<constraint_index>(-1)),
+        m_upper_bound_witness(static_cast<constraint_index>(-1)),
+        m_i(ri) {}
     ul_pair(const ul_pair & o): m_low_bound_witness(o.m_low_bound_witness), m_upper_bound_witness(o.m_upper_bound_witness), m_i(o.m_i) {}
 };