diff --git a/src/smt/theory_lra.cpp b/src/smt/theory_lra.cpp new file mode 100644 index 000000000..745aa5da5 --- /dev/null +++ b/src/smt/theory_lra.cpp @@ -0,0 +1,2624 @@ +/*++ +Copyright (c) 2016 Microsoft Corporation + +Module Name: + + theory_lra.cpp + +Abstract: + + + +Author: + + Lev Nachmanson (levnach) 2016-25-3 + Nikolaj Bjorner (nbjorner) + +Revision History: + + +--*/ +#include "util/stopwatch.h" +#include "util/lp/lp_solver.h" +#include "util/lp/lp_primal_simplex.h" +#include "util/lp/lp_dual_simplex.h" +#include "util/lp/indexed_value.h" +#include "util/lp/lar_solver.h" +#include "util/nat_set.h" +#include "util/optional.h" +#include "lp_params.hpp" +#include "util/inf_rational.h" +#include "smt/smt_theory.h" +#include "smt/smt_context.h" +#include "smt/theory_lra.h" +#include "smt/proto_model/numeral_factory.h" +#include "smt/smt_model_generator.h" +#include "smt/arith_eq_adapter.h" +#include "util/nat_set.h" +#include "tactic/filter_model_converter.h" + +namespace lp { + enum bound_kind { lower_t, upper_t }; + + std::ostream& operator<<(std::ostream& out, bound_kind const& k) { + switch (k) { + case lower_t: return out << "<="; + case upper_t: return out << ">="; + } + return out; + } + + class bound { + smt::bool_var m_bv; + smt::theory_var m_var; + rational m_value; + bound_kind m_bound_kind; + + public: + bound(smt::bool_var bv, smt::theory_var v, rational const & val, bound_kind k): + m_bv(bv), + m_var(v), + m_value(val), + m_bound_kind(k) { + } + virtual ~bound() {} + smt::theory_var get_var() const { return m_var; } + smt::bool_var get_bv() const { return m_bv; } + bound_kind get_bound_kind() const { return m_bound_kind; } + rational const& get_value() const { return m_value; } + inf_rational get_value(bool is_true) const { + if (is_true) return inf_rational(m_value); // v >= value or v <= value + if (m_bound_kind == lower_t) return inf_rational(m_value, false); // v <= value - epsilon + return inf_rational(m_value, true); // v >= value + epsilon + } + virtual std::ostream& display(std::ostream& out) const { + return out << "v" << get_var() << " " << get_bound_kind() << " " << m_value; + } + }; + + std::ostream& operator<<(std::ostream& out, bound const& b) { + return b.display(out); + } + + struct stats { + unsigned m_assert_lower; + unsigned m_assert_upper; + unsigned m_add_rows; + unsigned m_bounds_propagations; + unsigned m_num_iterations; + unsigned m_num_iterations_with_no_progress; + unsigned m_num_factorizations; + unsigned m_need_to_solve_inf; + unsigned m_fixed_eqs; + unsigned m_conflicts; + unsigned m_bound_propagations1; + unsigned m_bound_propagations2; + unsigned m_assert_diseq; + unsigned m_make_feasible; + unsigned m_max_cols; + unsigned m_max_rows; + stats() { reset(); } + void reset() { + memset(this, 0, sizeof(*this)); + } + }; + + typedef optional opt_inf_rational; + + +} + +namespace smt { + + typedef ptr_vector lp_bounds; + + class theory_lra::imp { + + struct scope { + unsigned m_bounds_lim; + unsigned m_asserted_qhead; + unsigned m_asserted_atoms_lim; + unsigned m_delayed_terms_lim; + unsigned m_delayed_equalities_lim; + unsigned m_delayed_defs_lim; + unsigned m_underspecified_lim; + unsigned m_var_trail_lim; + expr* m_not_handled; + }; + + struct delayed_atom { + unsigned m_bv; + bool m_is_true; + delayed_atom(unsigned b, bool t): m_bv(b), m_is_true(t) {} + }; + + class resource_limit : public lean::lp_resource_limit { + imp& m_imp; + public: + resource_limit(imp& i): m_imp(i) { } + virtual bool get_cancel_flag() { return m_imp.m.canceled(); } + }; + + + 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 m_columns; + int m_print_counter = 0; + + // temporary values kept during internalization + struct internalize_state { + expr_ref_vector m_terms; + vector m_coeffs; + svector m_vars; + rational m_coeff; + ptr_vector m_terms_to_internalize; + internalize_state(ast_manager& m): m_terms(m) {} + void reset() { + m_terms.reset(); + m_coeffs.reset(); + m_coeff.reset(); + m_vars.reset(); + m_terms_to_internalize.reset(); + } + }; + ptr_vector m_internalize_states; + unsigned m_internalize_head; + + class scoped_internalize_state { + imp& m_imp; + internalize_state& m_st; + + internalize_state& push_internalize(imp& i) { + if (i.m_internalize_head == i.m_internalize_states.size()) { + i.m_internalize_states.push_back(alloc(internalize_state, i.m)); + } + internalize_state& st = *i.m_internalize_states[i.m_internalize_head++]; + st.reset(); + return st; + } + public: + scoped_internalize_state(imp& i): m_imp(i), m_st(push_internalize(i)) {} + ~scoped_internalize_state() { --m_imp.m_internalize_head; } + expr_ref_vector& terms() { return m_st.m_terms; } + vector& coeffs() { return m_st.m_coeffs; } + svector& vars() { return m_st.m_vars; } + rational& coeff() { return m_st.m_coeff; } + ptr_vector& terms_to_internalize() { return m_st.m_terms_to_internalize; } + void push(expr* e, rational c) { m_st.m_terms.push_back(e); m_st.m_coeffs.push_back(c); } + void set_back(unsigned i) { + if (terms().size() == i + 1) return; + terms()[i] = terms().back(); + coeffs()[i] = coeffs().back(); + terms().pop_back(); + coeffs().pop_back(); + } + }; + + typedef vector> var_coeffs; + struct delayed_def { + vector m_coeffs; + svector m_vars; + rational m_coeff; + theory_var m_var; + delayed_def(svector const& vars, vector const& coeffs, rational const& r, theory_var v): + m_coeffs(coeffs), m_vars(vars), m_coeff(r), m_var(v) {} + }; + + svector m_theory_var2var_index; // translate from theory variables to lar vars + svector m_var_index2theory_var; // reverse map from lp_solver variables to theory variables + svector m_term_index2theory_var; // reverse map from lp_solver variables to theory variables + var_coeffs m_left_side; // constraint left side + mutable std::unordered_map m_variable_values; // current model + + enum constraint_source { + inequality_source, + equality_source, + definition_source, + null_source + }; + svector m_constraint_sources; + svector m_inequalities; // asserted rows corresponding to inequality literals. + svector m_equalities; // asserted rows corresponding to equalities. + svector m_definitions; // asserted rows corresponding to definitions + + bool m_delay_constraints; // configuration + svector m_asserted_atoms; + app_ref_vector m_delayed_terms; + svector> m_delayed_equalities; + vector m_delayed_defs; + expr* m_not_handled; + ptr_vector m_underspecified; + unsigned_vector m_var_trail; + vector > m_use_list; // bounds where variables are used. + + // attributes for incremental version: + u_map m_bool_var2bound; + vector m_bounds; + unsigned_vector m_unassigned_bounds; + unsigned_vector m_bounds_trail; + unsigned m_asserted_qhead; + + svector m_to_check; // rows that should be checked for theory propagation + + svector > m_assume_eq_candidates; + unsigned m_assume_eq_head; + + unsigned m_num_conflicts; + + + struct var_value_eq { + imp & m_th; + var_value_eq(imp & th):m_th(th) {} + bool operator()(theory_var v1, theory_var v2) const { return m_th.get_ivalue(v1) == m_th.get_ivalue(v2) && m_th.is_int(v1) == m_th.is_int(v2); } + }; + struct var_value_hash { + imp & m_th; + var_value_hash(imp & th):m_th(th) {} + unsigned operator()(theory_var v) const { return std::hash()(m_th.get_ivalue(v)); } + }; + int_hashtable m_model_eqs; + + + svector m_scopes; + lp::stats m_stats; + arith_factory* m_factory; + scoped_ptr m_solver; + resource_limit m_resource_limit; + lp_bounds m_new_bounds; + + + context& ctx() const { return th.get_context(); } + theory_id get_id() const { return th.get_id(); } + bool is_int(theory_var v) const { return is_int(get_enode(v)); } + bool is_int(enode* n) const { return a.is_int(n->get_owner()); } + enode* get_enode(theory_var v) const { return th.get_enode(v); } + enode* get_enode(expr* e) const { return ctx().get_enode(e); } + expr* get_owner(theory_var v) const { return get_enode(v)->get_owner(); } + + void init_solver() { + 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(m_lp_params.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->settings().set_ostream(0); + } + + void found_not_handled(expr* n) { + m_not_handled = n; + if (is_app(n) && is_underspecified(to_app(n))) { + m_underspecified.push_back(to_app(n)); + } + TRACE("arith", tout << "Unhandled: " << mk_pp(n, m) << "\n";); + } + + bool is_numeral(expr* term, rational& r) { + rational mul(1); + do { + if (a.is_numeral(term, r)) { + r *= mul; + return true; + } + if (a.is_uminus(term, term)) { + mul.neg(); + continue; + } + if (a.is_to_real(term, term)) { + continue; + } + return false; + } + while (false); + return false; + } + + void linearize_term(expr* term, scoped_internalize_state& st) { + st.push(term, rational::one()); + linearize(st); + } + + void linearize_ineq(expr* lhs, expr* rhs, scoped_internalize_state& st) { + st.push(lhs, rational::one()); + st.push(rhs, rational::minus_one()); + linearize(st); + } + + void linearize(scoped_internalize_state& st) { + expr_ref_vector & terms = st.terms(); + svector& vars = st.vars(); + vector& coeffs = st.coeffs(); + rational& coeff = st.coeff(); + rational r; + expr* n1, *n2; + unsigned index = 0; + while (index < terms.size()) { + SASSERT(index >= vars.size()); + expr* n = terms[index].get(); + st.terms_to_internalize().push_back(n); + if (a.is_add(n)) { + unsigned sz = to_app(n)->get_num_args(); + for (unsigned i = 0; i < sz; ++i) { + st.push(to_app(n)->get_arg(i), coeffs[index]); + } + st.set_back(index); + } + else if (a.is_sub(n)) { + unsigned sz = to_app(n)->get_num_args(); + terms[index] = to_app(n)->get_arg(0); + for (unsigned i = 1; i < sz; ++i) { + st.push(to_app(n)->get_arg(i), -coeffs[index]); + } + } + else if (a.is_mul(n, n1, n2) && is_numeral(n1, r)) { + coeffs[index] *= r; + terms[index] = n2; + st.terms_to_internalize().push_back(n1); + } + else if (a.is_mul(n, n1, n2) && is_numeral(n2, r)) { + coeffs[index] *= r; + terms[index] = n1; + st.terms_to_internalize().push_back(n2); + } + else if (a.is_numeral(n, r)) { + coeff += coeffs[index]*r; + ++index; + } + else if (a.is_uminus(n, n1)) { + coeffs[index].neg(); + terms[index] = n1; + } + else if (is_app(n) && a.get_family_id() == to_app(n)->get_family_id()) { + app* t = to_app(n); + found_not_handled(n); + internalize_args(t); + mk_enode(t); + theory_var v = mk_var(n); + coeffs[vars.size()] = coeffs[index]; + vars.push_back(v); + ++index; + } + else { + if (is_app(n)) { + internalize_args(to_app(n)); + } + if (a.is_int(n)) { + found_not_handled(n); + } + theory_var v = mk_var(n); + coeffs[vars.size()] = coeffs[index]; + vars.push_back(v); + ++index; + } + } + for (unsigned i = st.terms_to_internalize().size(); i > 0; ) { + --i; + expr* n = st.terms_to_internalize()[i]; + if (is_app(n)) { + mk_enode(to_app(n)); + } + } + st.terms_to_internalize().reset(); + } + + void internalize_args(app* t) { + for (unsigned i = 0; reflect(t) && i < t->get_num_args(); ++i) { + if (!ctx().e_internalized(t->get_arg(i))) { + ctx().internalize(t->get_arg(i), false); + } + } + } + + enode * mk_enode(app * n) { + if (ctx().e_internalized(n)) { + return get_enode(n); + } + else { + return ctx().mk_enode(n, !reflect(n), false, enable_cgc_for(n)); + } + } + + bool enable_cgc_for(app * n) const { + // Congruence closure is not enabled for (+ ...) and (* ...) applications. + return !(n->get_family_id() == get_id() && (n->get_decl_kind() == OP_ADD || n->get_decl_kind() == OP_MUL)); + } + + + void mk_clause(literal l1, literal l2, unsigned num_params, parameter * params) { + TRACE("arith", literal lits[2]; lits[0] = l1; lits[1] = l2; ctx().display_literals_verbose(tout, 2, lits); tout << "\n";); + ctx().mk_th_axiom(get_id(), l1, l2, num_params, params); + } + + void mk_clause(literal l1, literal l2, literal l3, unsigned num_params, parameter * params) { + TRACE("arith", literal lits[3]; lits[0] = l1; lits[1] = l2; lits[2] = l3; ctx().display_literals_verbose(tout, 3, lits); tout << "\n";); + ctx().mk_th_axiom(get_id(), l1, l2, l3, num_params, params); + } + + bool is_underspecified(app* n) const { + if (n->get_family_id() == get_id()) { + switch (n->get_decl_kind()) { + case OP_DIV: + case OP_IDIV: + case OP_REM: + case OP_MOD: + return true; + default: + break; + } + } + return false; + } + + bool reflect(app* n) const { + return m_arith_params.m_arith_reflect || is_underspecified(n); + } + + theory_var mk_var(expr* n, bool internalize = true) { + if (!ctx().e_internalized(n)) { + ctx().internalize(n, false); + } + enode* e = get_enode(n); + theory_var v; + if (!th.is_attached_to_var(e)) { + v = th.mk_var(e); + SASSERT(m_bounds.size() <= static_cast(v) || m_bounds[v].empty()); + if (m_bounds.size() <= static_cast(v)) { + m_bounds.push_back(lp_bounds()); + m_unassigned_bounds.push_back(0); + } + ctx().attach_th_var(e, &th, v); + } + else { + v = e->get_th_var(get_id()); + } + SASSERT(null_theory_var != v); + return v; + } + + lean::var_index get_var_index(theory_var v) { + lean::var_index result = UINT_MAX; + if (m_theory_var2var_index.size() > static_cast(v)) { + result = m_theory_var2var_index[v]; + } + if (result == UINT_MAX) { + result = m_solver->add_var(v); + m_theory_var2var_index.setx(v, result, UINT_MAX); + m_var_index2theory_var.setx(result, v, UINT_MAX); + m_var_trail.push_back(v); + } + return result; + } + + void init_left_side(scoped_internalize_state& st) { + SASSERT(all_zeros(m_columns)); + svector const& vars = st.vars(); + vector const& coeffs = st.coeffs(); + for (unsigned i = 0; i < vars.size(); ++i) { + theory_var var = vars[i]; + rational const& coeff = coeffs[i]; + if (m_columns.size() <= static_cast(var)) { + m_columns.setx(var, coeff, rational::zero()); + } + else { + m_columns[var] += coeff; + } + } + m_left_side.clear(); + // reset the coefficients after they have been used. + for (unsigned i = 0; i < vars.size(); ++i) { + theory_var var = vars[i]; + rational const& r = m_columns[var]; + if (!r.is_zero()) { + m_left_side.push_back(std::make_pair(r, get_var_index(var))); + m_columns[var].reset(); + } + } + SASSERT(all_zeros(m_columns)); + } + + bool all_zeros(vector const& v) const { + for (unsigned i = 0; i < v.size(); ++i) { + if (!v[i].is_zero()) { + return false; + } + } + return true; + } + + void add_eq_constraint(lean::constraint_index index, enode* n1, enode* n2) { + m_constraint_sources.setx(index, equality_source, null_source); + m_equalities.setx(index, enode_pair(n1, n2), enode_pair(0, 0)); + ++m_stats.m_add_rows; + } + + void add_ineq_constraint(lean::constraint_index index, literal lit) { + m_constraint_sources.setx(index, inequality_source, null_source); + m_inequalities.setx(index, lit, null_literal); + ++m_stats.m_add_rows; + TRACE("arith", m_solver->print_constraint(index, tout); tout << "\n";); + } + + void add_def_constraint(lean::constraint_index index, theory_var v) { + m_constraint_sources.setx(index, definition_source, null_source); + m_definitions.setx(index, v, null_theory_var); + ++m_stats.m_add_rows; + } + + void internalize_eq(delayed_def const& d) { + scoped_internalize_state st(*this); + st.vars().append(d.m_vars); + st.coeffs().append(d.m_coeffs); + init_left_side(st); + add_def_constraint(m_solver->add_constraint(m_left_side, lean::EQ, -d.m_coeff), d.m_var); + } + + void internalize_eq(theory_var v1, theory_var v2) { + enode* n1 = get_enode(v1); + enode* n2 = get_enode(v2); + scoped_internalize_state st(*this); + st.vars().push_back(v1); + st.vars().push_back(v2); + st.coeffs().push_back(rational::one()); + st.coeffs().push_back(rational::minus_one()); + init_left_side(st); + add_eq_constraint(m_solver->add_constraint(m_left_side, lean::EQ, rational::zero()), n1, n2); + TRACE("arith", + tout << "v" << v1 << " = " << "v" << v2 << ": " + << mk_pp(n1->get_owner(), m) << " = " << mk_pp(n2->get_owner(), m) << "\n";); + } + + void del_bounds(unsigned old_size) { + for (unsigned i = m_bounds_trail.size(); i > old_size; ) { + --i; + unsigned v = m_bounds_trail[i]; + lp::bound* b = m_bounds[v].back(); + // del_use_lists(b); + dealloc(b); + m_bounds[v].pop_back(); + } + m_bounds_trail.shrink(old_size); + } + + void updt_unassigned_bounds(theory_var v, int inc) { + TRACE("arith", tout << "v" << v << " " << m_unassigned_bounds[v] << " += " << inc << "\n";); + ctx().push_trail(vector_value_trail(m_unassigned_bounds, v)); + m_unassigned_bounds[v] += inc; + } + + bool is_unit_var(scoped_internalize_state& st) { + return st.coeff().is_zero() && st.vars().size() == 1 && st.coeffs()[0].is_one(); + } + + theory_var internalize_def(app* term, scoped_internalize_state& st) { + linearize_term(term, st); + if (is_unit_var(st)) { + return st.vars()[0]; + } + else { + theory_var v = mk_var(term); + SASSERT(null_theory_var != v); + st.coeffs().resize(st.vars().size() + 1); + st.coeffs()[st.vars().size()] = rational::minus_one(); + st.vars().push_back(v); + return v; + } + } + + // term - v = 0 + theory_var internalize_def(app* term) { + scoped_internalize_state st(*this); + linearize_term(term, st); + if (is_unit_var(st)) { + return st.vars()[0]; + } + else { + init_left_side(st); + theory_var v = mk_var(term); + lean::var_index vi = m_theory_var2var_index.get(v, UINT_MAX); + if (vi == UINT_MAX) { + vi = m_solver->add_term(m_left_side, st.coeff()); + m_theory_var2var_index.setx(v, vi, UINT_MAX); + if (m_solver->is_term(vi)) { + m_term_index2theory_var.setx(m_solver->adjust_term_index(vi), v, UINT_MAX); + } + else { + m_var_index2theory_var.setx(vi, v, UINT_MAX); + } + m_var_trail.push_back(v); + TRACE("arith", tout << "v" << v << " := " << mk_pp(term, m) << " slack: " << vi << " scopes: " << m_scopes.size() << "\n"; + m_solver->print_term(m_solver->get_term(vi), tout); tout << "\n";); + } + rational val; + if (a.is_numeral(term, val)) { + m_fixed_var_table.insert(value_sort_pair(val, is_int(v)), v); + } + return v; + } + } + + + public: + imp(theory_lra& th, ast_manager& m, theory_arith_params& p): + th(th), m(m), m_arith_params(p), a(m), + m_arith_eq_adapter(th, p, a), + m_internalize_head(0), + m_delay_constraints(false), + m_delayed_terms(m), + m_not_handled(0), + m_asserted_qhead(0), + 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_resource_limit(*this) { + init_solver(); + } + + ~imp() { + del_bounds(0); + std::for_each(m_internalize_states.begin(), m_internalize_states.end(), delete_proc()); + } + + bool internalize_atom(app * atom, bool gate_ctx) { + if (m_delay_constraints) { + return internalize_atom_lazy(atom, gate_ctx); + } + else { + return internalize_atom_strict(atom, gate_ctx); + } + } + + bool internalize_atom_strict(app * atom, bool gate_ctx) { + SASSERT(!ctx().b_internalized(atom)); + bool_var bv = ctx().mk_bool_var(atom); + ctx().set_var_theory(bv, get_id()); + expr* n1, *n2; + rational r; + lp::bound_kind k; + theory_var v = null_theory_var; + if (a.is_le(atom, n1, n2) && is_numeral(n2, r) && is_app(n1)) { + v = internalize_def(to_app(n1)); + k = lp::upper_t; + } + else if (a.is_ge(atom, n1, n2) && is_numeral(n2, r) && is_app(n1)) { + v = internalize_def(to_app(n1)); + k = lp::lower_t; + } + else { + TRACE("arith", tout << "Could not internalize " << mk_pp(atom, m) << "\n";); + found_not_handled(atom); + return true; + } + lp::bound* b = alloc(lp::bound, bv, v, r, k); + m_bounds[v].push_back(b); + updt_unassigned_bounds(v, +1); + m_bounds_trail.push_back(v); + m_bool_var2bound.insert(bv, b); + TRACE("arith", tout << "Internalized " << mk_pp(atom, m) << "\n";); + mk_bound_axioms(*b); + //add_use_lists(b); + return true; + } + + bool internalize_atom_lazy(app * atom, bool gate_ctx) { + SASSERT(!ctx().b_internalized(atom)); + bool_var bv = ctx().mk_bool_var(atom); + ctx().set_var_theory(bv, get_id()); + expr* n1, *n2; + rational r; + lp::bound_kind k; + theory_var v = null_theory_var; + scoped_internalize_state st(*this); + if (a.is_le(atom, n1, n2) && is_numeral(n2, r) && is_app(n1)) { + v = internalize_def(to_app(n1), st); + k = lp::upper_t; + } + else if (a.is_ge(atom, n1, n2) && is_numeral(n2, r) && is_app(n1)) { + v = internalize_def(to_app(n1), st); + k = lp::lower_t; + } + else { + TRACE("arith", tout << "Could not internalize " << mk_pp(atom, m) << "\n";); + found_not_handled(atom); + return true; + } + lp::bound* b = alloc(lp::bound, bv, v, r, k); + m_bounds[v].push_back(b); + updt_unassigned_bounds(v, +1); + m_bounds_trail.push_back(v); + m_bool_var2bound.insert(bv, b); + TRACE("arith", tout << "Internalized " << mk_pp(atom, m) << "\n";); + if (!is_unit_var(st) && m_bounds[v].size() == 1) { + m_delayed_defs.push_back(delayed_def(st.vars(), st.coeffs(), st.coeff(), v)); + } + return true; + } + + bool internalize_term(app * term) { + if (ctx().e_internalized(term) && th.is_attached_to_var(ctx().get_enode(term))) { + // skip + } + else if (m_delay_constraints) { + scoped_internalize_state st(*this); + linearize_term(term, st); // ensure that a theory_var was created. + SASSERT(ctx().e_internalized(term)); + if(!th.is_attached_to_var(ctx().get_enode(term))) { + mk_var(term); + } + m_delayed_terms.push_back(term); + } + else { + internalize_def(term); + } + return true; + } + + void internalize_eq_eh(app * atom, bool_var) { + expr* lhs, *rhs; + VERIFY(m.is_eq(atom, lhs, rhs)); + enode * n1 = get_enode(lhs); + enode * n2 = get_enode(rhs); + if (n1->get_th_var(get_id()) != null_theory_var && + n2->get_th_var(get_id()) != null_theory_var && + n1 != n2) { + TRACE("arith", tout << mk_pp(atom, m) << "\n";); + m_arith_eq_adapter.mk_axioms(n1, n2); + } + } + + void assign_eh(bool_var v, bool is_true) { + TRACE("arith", tout << mk_pp(ctx().bool_var2expr(v), m) << " " << (is_true?"true":"false") << "\n";); + m_asserted_atoms.push_back(delayed_atom(v, is_true)); + } + + void new_eq_eh(theory_var v1, theory_var v2) { + if (m_delay_constraints) { + m_delayed_equalities.push_back(std::make_pair(v1, v2)); + } + else { + // or internalize_eq(v1, v2); + m_arith_eq_adapter.new_eq_eh(v1, v2); + } + } + + bool use_diseqs() const { + return true; + } + + void new_diseq_eh(theory_var v1, theory_var v2) { + TRACE("arith", tout << "v" << v1 << " != " << "v" << v2 << "\n";); + ++m_stats.m_assert_diseq; + m_arith_eq_adapter.new_diseq_eh(v1, v2); + } + + void push_scope_eh() { + m_scopes.push_back(scope()); + scope& s = m_scopes.back(); + s.m_bounds_lim = m_bounds_trail.size(); + s.m_asserted_qhead = m_asserted_qhead; + s.m_asserted_atoms_lim = m_asserted_atoms.size(); + s.m_delayed_terms_lim = m_delayed_terms.size(); + s.m_delayed_equalities_lim = m_delayed_equalities.size(); + s.m_delayed_defs_lim = m_delayed_defs.size(); + s.m_not_handled = m_not_handled; + s.m_underspecified_lim = m_underspecified.size(); + s.m_var_trail_lim = m_var_trail.size(); + if (!m_delay_constraints) m_solver->push(); + } + + void pop_scope_eh(unsigned num_scopes) { + if (num_scopes == 0) { + return; + } + unsigned old_size = m_scopes.size() - num_scopes; + del_bounds(m_scopes[old_size].m_bounds_lim); + for (unsigned i = m_scopes[old_size].m_var_trail_lim; i < m_var_trail.size(); ++i) { + lean::var_index vi = m_theory_var2var_index[m_var_trail[i]]; + if (m_solver->is_term(vi)) { + unsigned ti = m_solver->adjust_term_index(vi); + m_term_index2theory_var[ti] = UINT_MAX; + } + else if (vi < m_var_index2theory_var.size()) { + m_var_index2theory_var[vi] = UINT_MAX; + } + m_theory_var2var_index[m_var_trail[i]] = UINT_MAX; + } + m_asserted_atoms.shrink(m_scopes[old_size].m_asserted_atoms_lim); + m_delayed_terms.shrink(m_scopes[old_size].m_delayed_terms_lim); + m_delayed_defs.shrink(m_scopes[old_size].m_delayed_defs_lim); + m_delayed_equalities.shrink(m_scopes[old_size].m_delayed_equalities_lim); + m_asserted_qhead = m_scopes[old_size].m_asserted_qhead; + m_underspecified.shrink(m_scopes[old_size].m_underspecified_lim); + m_var_trail.shrink(m_scopes[old_size].m_var_trail_lim); + m_not_handled = m_scopes[old_size].m_not_handled; + m_scopes.resize(old_size); + if (!m_delay_constraints) m_solver->pop(num_scopes); + // VERIFY(l_false != make_feasible()); + m_new_bounds.reset(); + m_to_check.reset(); + TRACE("arith", tout << "num scopes: " << num_scopes << " new scope level: " << m_scopes.size() << "\n";); + } + + void restart_eh() { + m_arith_eq_adapter.restart_eh(); + } + + void relevant_eh(app* n) { + TRACE("arith", tout << mk_pp(n, m) << "\n";); + expr* n1, *n2; + if (a.is_mod(n, n1, n2)) + mk_idiv_mod_axioms(n1, n2); + else if (a.is_rem(n, n1, n2)) + mk_rem_axiom(n1, n2); + else if (a.is_div(n, n1, n2)) + mk_div_axiom(n1, n2); + else if (a.is_to_int(n)) + mk_to_int_axiom(n); + else if (a.is_is_int(n)) + mk_is_int_axiom(n); + } + + // n < 0 || rem(a, n) = mod(a, n) + // !n < 0 || rem(a, n) = -mod(a, n) + void mk_rem_axiom(expr* dividend, expr* divisor) { + expr_ref zero(a.mk_int(0), m); + expr_ref rem(a.mk_rem(dividend, divisor), m); + expr_ref mod(a.mk_mod(dividend, divisor), m); + expr_ref mmod(a.mk_uminus(mod), m); + literal dgez = mk_literal(a.mk_ge(divisor, zero)); + mk_axiom(~dgez, th.mk_eq(rem, mod, false)); + mk_axiom( dgez, th.mk_eq(rem, mmod, false)); + } + + // q = 0 or q * (p div q) = p + void mk_div_axiom(expr* p, expr* q) { + if (a.is_zero(q)) return; + literal eqz = th.mk_eq(q, a.mk_real(0), false); + literal eq = th.mk_eq(a.mk_mul(q, a.mk_div(p, q)), p, false); + mk_axiom(eqz, eq); + } + + // to_int (to_real x) = x + // to_real(to_int(x)) <= x < to_real(to_int(x)) + 1 + void mk_to_int_axiom(app* n) { + expr* x, *y; + VERIFY (a.is_to_int(n, x)); + if (a.is_to_real(x, y)) { + mk_axiom(th.mk_eq(y, n, false)); + } + else { + expr_ref to_r(a.mk_to_real(n), m); + expr_ref lo(a.mk_le(a.mk_sub(to_r, x), a.mk_real(0)), m); + expr_ref hi(a.mk_ge(a.mk_sub(x, to_r), a.mk_real(1)), m); + mk_axiom(mk_literal(lo)); + mk_axiom(~mk_literal(hi)); + } + } + + // is_int(x) <=> to_real(to_int(x)) = x + void mk_is_int_axiom(app* n) { + expr* x; + VERIFY(a.is_is_int(n, x)); + literal eq = th.mk_eq(a.mk_to_real(a.mk_to_int(x)), x, false); + literal is_int = ctx().get_literal(n); + mk_axiom(~is_int, eq); + mk_axiom(is_int, ~eq); + } + + void mk_idiv_mod_axioms(expr * p, expr * q) { + if (a.is_zero(q)) { + return; + } + // if q is zero, then idiv and mod are uninterpreted functions. + expr_ref div(a.mk_idiv(p, q), m); + expr_ref mod(a.mk_mod(p, q), m); + expr_ref zero(a.mk_int(0), m); + literal q_ge_0 = mk_literal(a.mk_ge(q, zero)); + literal q_le_0 = mk_literal(a.mk_le(q, zero)); + // literal eqz = th.mk_eq(q, zero, false); + literal eq = th.mk_eq(a.mk_add(a.mk_mul(q, div), mod), p, false); + literal mod_ge_0 = mk_literal(a.mk_ge(mod, zero)); + // q >= 0 or p = (p mod q) + q * (p div q) + // q <= 0 or p = (p mod q) + q * (p div q) + // q >= 0 or (p mod q) >= 0 + // q <= 0 or (p mod q) >= 0 + // q <= 0 or (p mod q) < q + // q >= 0 or (p mod q) < -q + mk_axiom(q_ge_0, eq); + mk_axiom(q_le_0, eq); + mk_axiom(q_ge_0, mod_ge_0); + mk_axiom(q_le_0, mod_ge_0); + mk_axiom(q_le_0, ~mk_literal(a.mk_ge(a.mk_sub(mod, q), zero))); + mk_axiom(q_ge_0, ~mk_literal(a.mk_ge(a.mk_add(mod, q), zero))); + rational k; + if (m_arith_params.m_arith_enum_const_mod && a.is_numeral(q, k) && + k.is_pos() && k < rational(8)) { + unsigned _k = k.get_unsigned(); + literal_buffer lits; + for (unsigned j = 0; j < _k; ++j) { + literal mod_j = th.mk_eq(mod, a.mk_int(j), false); + lits.push_back(mod_j); + ctx().mark_as_relevant(mod_j); + } + ctx().mk_th_axiom(get_id(), lits.size(), lits.begin()); + } + } + + void mk_axiom(literal l) { + ctx().mk_th_axiom(get_id(), false_literal, l); + if (ctx().relevancy()) { + ctx().mark_as_relevant(l); + } + } + + void mk_axiom(literal l1, literal l2) { + if (l1 == false_literal) { + mk_axiom(l2); + return; + } + ctx().mk_th_axiom(get_id(), l1, l2); + if (ctx().relevancy()) { + ctx().mark_as_relevant(l1); + expr_ref e(m); + ctx().literal2expr(l2, e); + ctx().add_rel_watch(~l1, e); + } + } + + void mk_axiom(literal l1, literal l2, literal l3) { + ctx().mk_th_axiom(get_id(), l1, l2, l3); + if (ctx().relevancy()) { + expr_ref e(m); + ctx().mark_as_relevant(l1); + ctx().literal2expr(l2, e); + ctx().add_rel_watch(~l1, e); + ctx().literal2expr(l3, e); + ctx().add_rel_watch(~l2, e); + } + } + + literal mk_literal(expr* e) { + expr_ref pinned(e, m); + if (!ctx().e_internalized(e)) { + ctx().internalize(e, false); + } + return ctx().get_literal(e); + } + + + void init_search_eh() { + m_arith_eq_adapter.init_search_eh(); + m_num_conflicts = 0; + } + + bool can_get_value(theory_var v) const { + return + (v != null_theory_var) && + (v < static_cast(m_theory_var2var_index.size())) && + (UINT_MAX != m_theory_var2var_index[v]) && + (m_solver->is_term(m_theory_var2var_index[v]) || m_variable_values.count(m_theory_var2var_index[v]) > 0); + } + + + bool can_get_ivalue(theory_var v) const { + if (v == null_theory_var || (v >= static_cast(m_theory_var2var_index.size()))) + return false; + return m_solver->var_is_registered(m_theory_var2var_index[v]); + } + + lean::impq get_ivalue(theory_var v) const { + lean_assert(can_get_ivalue(v)); + lean::var_index vi = m_theory_var2var_index[v]; + if (!m_solver->is_term(vi)) + return m_solver->get_value(vi); + + const lean::lar_term& term = m_solver->get_term(vi); + lean::impq result(term.m_v); + for (const auto & i: term.m_coeffs) { + result += m_solver->get_value(i.first) * i.second; + } + return result; + } + + + rational get_value(theory_var v) const { + if (!can_get_value(v)) return rational::zero(); + lean::var_index vi = m_theory_var2var_index[v]; + if (m_variable_values.count(vi) > 0) { + return m_variable_values[vi]; + } + if (m_solver->is_term(vi)) { + const lean::lar_term& term = m_solver->get_term(vi); + rational result = term.m_v; + for (auto i = term.m_coeffs.begin(); i != term.m_coeffs.end(); ++i) { + result += m_variable_values[i->first] * i->second; + } + m_variable_values[vi] = result; + return result; + } + UNREACHABLE(); + return m_variable_values[vi]; + } + + void init_variable_values() { + if (m_solver.get() && th.get_num_vars() > 0) { + m_solver->get_model(m_variable_values); + } + } + + void reset_variable_values() { + m_variable_values.clear(); + } + + bool assume_eqs() { + svector vars; + theory_var sz = static_cast(th.get_num_vars()); + for (theory_var v = 0; v < sz; ++v) { + if (th.is_relevant_and_shared(get_enode(v))) { + vars.push_back(m_theory_var2var_index[v]); + } + } + if (vars.empty()) { + return false; + } + TRACE("arith", + for (theory_var v = 0; v < sz; ++v) { + if (th.is_relevant_and_shared(get_enode(v))) { + tout << "v" << v << " " << m_theory_var2var_index[v] << " "; + } + } + tout << "\n"; + ); + m_solver->random_update(vars.size(), vars.c_ptr()); + m_model_eqs.reset(); + TRACE("arith", display(tout);); + + unsigned old_sz = m_assume_eq_candidates.size(); + bool result = false; + int start = ctx().get_random_value(); + for (theory_var i = 0; i < sz; ++i) { + theory_var v = (i + start) % sz; + enode* n1 = get_enode(v); + if (!th.is_relevant_and_shared(n1)) { + continue; + } + if (!can_get_ivalue(v)) { + continue; + } + theory_var other = m_model_eqs.insert_if_not_there(v); + if (other == v) { + continue; + } + enode* n2 = get_enode(other); + if (n1->get_root() != n2->get_root()) { + TRACE("arith", tout << mk_pp(n1->get_owner(), m) << " = " << mk_pp(n2->get_owner(), m) << "\n"; + tout << mk_pp(n1->get_owner(), m) << " = " << mk_pp(n2->get_owner(), m) << "\n"; + tout << "v" << v << " = " << "v" << other << "\n";); + m_assume_eq_candidates.push_back(std::make_pair(v, other)); + result = true; + } + } + + if (result) { + ctx().push_trail(restore_size_trail, false>(m_assume_eq_candidates, old_sz)); + } + + return delayed_assume_eqs(); + } + + bool delayed_assume_eqs() { + if (m_assume_eq_head == m_assume_eq_candidates.size()) + return false; + + ctx().push_trail(value_trail(m_assume_eq_head)); + while (m_assume_eq_head < m_assume_eq_candidates.size()) { + std::pair const & p = m_assume_eq_candidates[m_assume_eq_head]; + theory_var v1 = p.first; + theory_var v2 = p.second; + enode* n1 = get_enode(v1); + enode* n2 = get_enode(v2); + m_assume_eq_head++; + CTRACE("arith", + get_ivalue(v1) == get_ivalue(v2) && n1->get_root() != n2->get_root(), + tout << "assuming eq: v" << v1 << " = v" << v2 << "\n";); + if (get_ivalue(v1) == get_ivalue(v2) && n1->get_root() != n2->get_root() && th.assume_eq(n1, n2)) { + return true; + } + } + return false; + } + + bool has_delayed_constraints() const { + return !(m_asserted_atoms.empty() && m_delayed_terms.empty() && m_delayed_equalities.empty()); + } + + final_check_status final_check_eh() { + lbool is_sat = l_true; + if (m_delay_constraints) { + init_solver(); + for (unsigned i = 0; i < m_asserted_atoms.size(); ++i) { + bool_var bv = m_asserted_atoms[i].m_bv; + assert_bound(bv, m_asserted_atoms[i].m_is_true, *m_bool_var2bound.find(bv)); + } + for (unsigned i = 0; i < m_delayed_terms.size(); ++i) { + internalize_def(m_delayed_terms[i].get()); + } + for (unsigned i = 0; i < m_delayed_defs.size(); ++i) { + internalize_eq(m_delayed_defs[i]); + } + for (unsigned i = 0; i < m_delayed_equalities.size(); ++i) { + std::pair const& eq = m_delayed_equalities[i]; + internalize_eq(eq.first, eq.second); + } + is_sat = make_feasible(); + } + else if (m_solver->get_status() != lean::lp_status::OPTIMAL) { + is_sat = make_feasible(); + } + switch (is_sat) { + case l_true: + if (delayed_assume_eqs()) { + return FC_CONTINUE; + } + if (assume_eqs()) { + return FC_CONTINUE; + } + if (m_not_handled != 0) { + return FC_GIVEUP; + } + return FC_DONE; + case l_false: + set_conflict(); + return FC_CONTINUE; + case l_undef: + return m.canceled() ? FC_CONTINUE : FC_GIVEUP; + default: + UNREACHABLE(); + break; + } + return FC_GIVEUP; + } + + + /** + \brief We must redefine this method, because theory of arithmetic contains + underspecified operators such as division by 0. + (/ a b) is essentially an uninterpreted function when b = 0. + Thus, 'a' must be considered a shared var if it is the child of an underspecified operator. + + if merge(a / b, x + y) and a / b is root, then x + y become shared and all z + u in equivalence class of x + y. + + + TBD: when the set of underspecified subterms is small, compute the shared variables below it. + Recompute the set if there are merges that invalidate it. + Use the set to determine if a variable is shared. + */ + bool is_shared(theory_var v) const { + if (m_underspecified.empty()) { + return false; + } + enode * n = get_enode(v); + enode * r = n->get_root(); + unsigned usz = m_underspecified.size(); + TRACE("shared", tout << ctx().get_scope_level() << " " << v << " " << r->get_num_parents() << "\n";); + if (r->get_num_parents() > 2*usz) { + for (unsigned i = 0; i < usz; ++i) { + app* u = m_underspecified[i]; + unsigned sz = u->get_num_args(); + for (unsigned j = 0; j < sz; ++j) { + if (ctx().get_enode(u->get_arg(j))->get_root() == r) { + return true; + } + } + } + } + else { + enode_vector::const_iterator it = r->begin_parents(); + enode_vector::const_iterator end = r->end_parents(); + for (; it != end; ++it) { + enode * parent = *it; + if (is_underspecified(parent->get_owner())) { + return true; + } + } + } + return false; + } + + bool can_propagate() { +#if 0 + if (ctx().at_base_level() && has_delayed_constraints()) { + // we could add the delayed constraints here directly to the tableau instead of using bounds variables. + } +#endif + return m_asserted_atoms.size() > m_asserted_qhead; + } + + void propagate() { + flush_bound_axioms(); + if (!can_propagate()) { + return; + } + unsigned qhead = m_asserted_qhead; + while (m_asserted_qhead < m_asserted_atoms.size() && !ctx().inconsistent()) { + bool_var bv = m_asserted_atoms[m_asserted_qhead].m_bv; + bool is_true = m_asserted_atoms[m_asserted_qhead].m_is_true; + +#if 1 + m_to_check.push_back(bv); +#else + propagate_bound(bv, is_true, b); +#endif + if (!m_delay_constraints) { + lp::bound& b = *m_bool_var2bound.find(bv); + assert_bound(bv, is_true, b); + } + + ++m_asserted_qhead; + } + if (m_delay_constraints || ctx().inconsistent()) { + m_to_check.reset(); + return; + } + /*for (; qhead < m_asserted_atoms.size() && !ctx().inconsistent(); ++qhead) { + bool_var bv = m_asserted_atoms[qhead].m_bv; + bool is_true = m_asserted_atoms[qhead].m_is_true; + lp::bound& b = *m_bool_var2bound.find(bv); + propagate_bound_compound(bv, is_true, b); + }*/ + + lbool lbl = make_feasible(); + + switch(lbl) { + case l_false: + TRACE("arith", tout << "propagation conflict\n";); + set_conflict(); + break; + case l_true: + propagate_basic_bounds(); + propagate_bounds_with_lp_solver(); + break; + case l_undef: + break; + } + + } + + void propagate_bounds_with_lp_solver() { + if (BP_NONE == propagation_mode()) { + return; + } + int num_of_p = m_solver->settings().st().m_num_of_implied_bounds; + local_bound_propagator bp(*this); + m_solver->propagate_bounds_for_touched_rows(bp); + if (m.canceled()) { + return; + } + int new_num_of_p = m_solver->settings().st().m_num_of_implied_bounds; + CTRACE("arith", new_num_of_p > num_of_p, tout << "found " << new_num_of_p << " implied bounds\n";); + if (m_solver->get_status() == lean::lp_status::INFEASIBLE) { + set_conflict(); + } + else { + for (size_t i = 0; !m.canceled() && !ctx().inconsistent() && i < bp.m_ibounds.size(); ++i) { + propagate_lp_solver_bound(bp.m_ibounds[i]); + } + } + } + + bool bound_is_interesting(unsigned vi, lean::lconstraint_kind kind, const rational & bval) const { + theory_var v; + if (m_solver->is_term(vi)) { + v = m_term_index2theory_var.get(m_solver->adjust_term_index(vi), null_theory_var); + } + else { + v = m_var_index2theory_var.get(vi, null_theory_var); + } + + if (v == null_theory_var) return false; + + if (m_unassigned_bounds[v] == 0 || m_bounds.size() <= static_cast(v)) { + TRACE("arith", tout << "return\n";); + return false; + } + lp_bounds const& bounds = m_bounds[v]; + for (unsigned i = 0; i < bounds.size(); ++i) { + lp::bound* b = bounds[i]; + if (ctx().get_assignment(b->get_bv()) != l_undef) { + continue; + } + literal lit = is_bound_implied(kind, bval, *b); + if (lit == null_literal) { + continue; + } + return true; + } + return false; + } + + struct local_bound_propagator: public lean::bound_propagator { + imp & m_imp; + local_bound_propagator(imp& i) : bound_propagator(*i.m_solver), m_imp(i) {} + + bool bound_is_interesting(unsigned j, lean::lconstraint_kind kind, const rational & v) { + return m_imp.bound_is_interesting(j, kind, v); + } + + virtual void consume(rational const& v, unsigned j) { + m_imp.set_evidence(j); + } + }; + + + void propagate_lp_solver_bound(lean::implied_bound& be) { + + theory_var v; + lean::var_index vi = be.m_j; + if (m_solver->is_term(vi)) { + v = m_term_index2theory_var.get(m_solver->adjust_term_index(vi), null_theory_var); + } + else { + v = m_var_index2theory_var.get(vi, null_theory_var); + } + + if (v == null_theory_var) return; + TRACE("arith", tout << "v" << v << " " << be.kind() << " " << be.m_bound << "\n"; + // if (m_unassigned_bounds[v] == 0) m_solver->print_bound_evidence(be, tout); + ); + + + if (m_unassigned_bounds[v] == 0 || m_bounds.size() <= static_cast(v)) { + TRACE("arith", tout << "return\n";); + return; + } + lp_bounds const& bounds = m_bounds[v]; + bool first = true; + for (unsigned i = 0; i < bounds.size(); ++i) { + lp::bound* b = bounds[i]; + if (ctx().get_assignment(b->get_bv()) != l_undef) { + continue; + } + literal lit = is_bound_implied(be.kind(), be.m_bound, *b); + if (lit == null_literal) { + continue; + } + + m_solver->settings().st().m_num_of_implied_bounds ++; + if (first) { + first = false; + m_core.reset(); + m_eqs.reset(); + m_params.reset(); + local_bound_propagator bp(*this); + m_solver->explain_implied_bound(be, bp); + } + CTRACE("arith", m_unassigned_bounds[v] == 0, tout << "missed bound\n";); + updt_unassigned_bounds(v, -1); + TRACE("arith", + ctx().display_literals_verbose(tout, m_core); + tout << "\n --> "; + ctx().display_literal_verbose(tout, lit); + tout << "\n"; + // display_evidence(tout, m_explanation); + m_solver->print_implied_bound(be, tout); + ); + DEBUG_CODE( + for (auto& lit : m_core) { + SASSERT(ctx().get_assignment(lit) == l_true); + }); + ++m_stats.m_bound_propagations1; + assign(lit); + } + } + + literal_vector m_core2; + + void assign(literal lit) { + SASSERT(validate_assign(lit)); + if (m_core.size() < small_lemma_size() && m_eqs.empty()) { + m_core2.reset(); + for (unsigned i = 0; i < m_core.size(); ++i) { + m_core2.push_back(~m_core[i]); + } + m_core2.push_back(lit); + justification * js = 0; + if (proofs_enabled()) { + js = alloc(theory_lemma_justification, get_id(), ctx(), m_core2.size(), m_core2.c_ptr(), + m_params.size(), m_params.c_ptr()); + } + ctx().mk_clause(m_core2.size(), m_core2.c_ptr(), js, CLS_AUX_LEMMA, 0); + } + else { + ctx().assign( + lit, ctx().mk_justification( + ext_theory_propagation_justification( + get_id(), ctx().get_region(), m_core.size(), m_core.c_ptr(), + m_eqs.size(), m_eqs.c_ptr(), lit, m_params.size(), m_params.c_ptr()))); + } + } + + literal is_bound_implied(lean::lconstraint_kind k, rational const& value, lp::bound const& b) const { + if ((k == lean::LE || k == lean::LT) && b.get_bound_kind() == lp::upper_t && value <= b.get_value()) { + // v <= value <= b.get_value() => v <= b.get_value() + return literal(b.get_bv(), false); + } + if ((k == lean::GE || k == lean::GT) && b.get_bound_kind() == lp::lower_t && b.get_value() <= value) { + // b.get_value() <= value <= v => b.get_value() <= v + return literal(b.get_bv(), false); + } + if (k == lean::LE && b.get_bound_kind() == lp::lower_t && value < b.get_value()) { + // v <= value < b.get_value() => v < b.get_value() + return literal(b.get_bv(), true); + } + if (k == lean::LT && b.get_bound_kind() == lp::lower_t && value <= b.get_value()) { + // v < value <= b.get_value() => v < b.get_value() + return literal(b.get_bv(), true); + } + if (k == lean::GE && b.get_bound_kind() == lp::upper_t && b.get_value() < value) { + // b.get_value() < value <= v => b.get_value() < v + return literal(b.get_bv(), true); + } + if (k == lean::GT && b.get_bound_kind() == lp::upper_t && b.get_value() <= value) { + // b.get_value() <= value < v => b.get_value() < v + return literal(b.get_bv(), true); + } + + return null_literal; + } + + void mk_bound_axioms(lp::bound& b) { + if (!ctx().is_searching()) { + // + // NB. We make an assumption that user push calls propagation + // before internal scopes are pushed. This flushes all newly + // asserted atoms into the right context. + // + m_new_bounds.push_back(&b); + return; + } + theory_var v = b.get_var(); + lp::bound_kind kind1 = b.get_bound_kind(); + rational const& k1 = b.get_value(); + lp_bounds & bounds = m_bounds[v]; + + lp::bound* end = 0; + lp::bound* lo_inf = end, *lo_sup = end; + lp::bound* hi_inf = end, *hi_sup = end; + + for (unsigned i = 0; i < bounds.size(); ++i) { + lp::bound& other = *bounds[i]; + if (&other == &b) continue; + if (b.get_bv() == other.get_bv()) continue; + lp::bound_kind kind2 = other.get_bound_kind(); + rational const& k2 = other.get_value(); + if (k1 == k2 && kind1 == kind2) { + // the bounds are equivalent. + continue; + } + + SASSERT(k1 != k2 || kind1 != kind2); + if (kind2 == lp::lower_t) { + if (k2 < k1) { + if (lo_inf == end || k2 > lo_inf->get_value()) { + lo_inf = &other; + } + } + else if (lo_sup == end || k2 < lo_sup->get_value()) { + lo_sup = &other; + } + } + else if (k2 < k1) { + if (hi_inf == end || k2 > hi_inf->get_value()) { + hi_inf = &other; + } + } + else if (hi_sup == end || k2 < hi_sup->get_value()) { + hi_sup = &other; + } + } + if (lo_inf != end) mk_bound_axiom(b, *lo_inf); + if (lo_sup != end) mk_bound_axiom(b, *lo_sup); + if (hi_inf != end) mk_bound_axiom(b, *hi_inf); + if (hi_sup != end) mk_bound_axiom(b, *hi_sup); + } + + + void mk_bound_axiom(lp::bound& b1, lp::bound& b2) { + theory_var v = b1.get_var(); + literal l1(b1.get_bv()); + literal l2(b2.get_bv()); + rational const& k1 = b1.get_value(); + rational const& k2 = b2.get_value(); + lp::bound_kind kind1 = b1.get_bound_kind(); + lp::bound_kind kind2 = b2.get_bound_kind(); + bool v_is_int = is_int(v); + SASSERT(v == b2.get_var()); + if (k1 == k2 && kind1 == kind2) return; + SASSERT(k1 != k2 || kind1 != kind2); + parameter coeffs[3] = { parameter(symbol("farkas")), + parameter(rational(1)), parameter(rational(1)) }; + + if (kind1 == lp::lower_t) { + if (kind2 == lp::lower_t) { + if (k2 <= k1) { + mk_clause(~l1, l2, 3, coeffs); + } + else { + mk_clause(l1, ~l2, 3, coeffs); + } + } + else if (k1 <= k2) { + // k1 <= k2, k1 <= x or x <= k2 + mk_clause(l1, l2, 3, coeffs); + } + else { + // k1 > hi_inf, k1 <= x => ~(x <= hi_inf) + mk_clause(~l1, ~l2, 3, coeffs); + if (v_is_int && k1 == k2 + rational(1)) { + // k1 <= x or x <= k1-1 + mk_clause(l1, l2, 3, coeffs); + } + } + } + else if (kind2 == lp::lower_t) { + if (k1 >= k2) { + // k1 >= lo_inf, k1 >= x or lo_inf <= x + mk_clause(l1, l2, 3, coeffs); + } + else { + // k1 < k2, k2 <= x => ~(x <= k1) + mk_clause(~l1, ~l2, 3, coeffs); + if (v_is_int && k1 == k2 - rational(1)) { + // x <= k1 or k1+l <= x + mk_clause(l1, l2, 3, coeffs); + } + + } + } + else { + // kind1 == A_UPPER, kind2 == A_UPPER + if (k1 >= k2) { + // k1 >= k2, x <= k2 => x <= k1 + mk_clause(l1, ~l2, 3, coeffs); + } + else { + // k1 <= hi_sup , x <= k1 => x <= hi_sup + mk_clause(~l1, l2, 3, coeffs); + } + } + } + + typedef lp_bounds::iterator iterator; + + void flush_bound_axioms() { + CTRACE("arith", !m_new_bounds.empty(), tout << "flush bound axioms\n";); + + while (!m_new_bounds.empty()) { + lp_bounds atoms; + atoms.push_back(m_new_bounds.back()); + m_new_bounds.pop_back(); + theory_var v = atoms.back()->get_var(); + for (unsigned i = 0; i < m_new_bounds.size(); ++i) { + if (m_new_bounds[i]->get_var() == v) { + atoms.push_back(m_new_bounds[i]); + m_new_bounds[i] = m_new_bounds.back(); + m_new_bounds.pop_back(); + --i; + } + } + CTRACE("arith_verbose", !atoms.empty(), + for (unsigned i = 0; i < atoms.size(); ++i) { + atoms[i]->display(tout); tout << "\n"; + }); + lp_bounds occs(m_bounds[v]); + + std::sort(atoms.begin(), atoms.end(), compare_bounds()); + std::sort(occs.begin(), occs.end(), compare_bounds()); + + iterator begin1 = occs.begin(); + iterator begin2 = occs.begin(); + iterator end = occs.end(); + begin1 = first(lp::lower_t, begin1, end); + begin2 = first(lp::upper_t, begin2, end); + + iterator lo_inf = begin1, lo_sup = begin1; + iterator hi_inf = begin2, hi_sup = begin2; + iterator lo_inf1 = begin1, lo_sup1 = begin1; + iterator hi_inf1 = begin2, hi_sup1 = begin2; + bool flo_inf, fhi_inf, flo_sup, fhi_sup; + ptr_addr_hashtable visited; + for (unsigned i = 0; i < atoms.size(); ++i) { + lp::bound* a1 = atoms[i]; + lo_inf1 = next_inf(a1, lp::lower_t, lo_inf, end, flo_inf); + hi_inf1 = next_inf(a1, lp::upper_t, hi_inf, end, fhi_inf); + lo_sup1 = next_sup(a1, lp::lower_t, lo_sup, end, flo_sup); + hi_sup1 = next_sup(a1, lp::upper_t, hi_sup, end, fhi_sup); + if (lo_inf1 != end) lo_inf = lo_inf1; + if (lo_sup1 != end) lo_sup = lo_sup1; + if (hi_inf1 != end) hi_inf = hi_inf1; + if (hi_sup1 != end) hi_sup = hi_sup1; + if (!flo_inf) lo_inf = end; + if (!fhi_inf) hi_inf = end; + if (!flo_sup) lo_sup = end; + if (!fhi_sup) hi_sup = end; + visited.insert(a1); + if (lo_inf1 != end && lo_inf != end && !visited.contains(*lo_inf)) mk_bound_axiom(*a1, **lo_inf); + if (lo_sup1 != end && lo_sup != end && !visited.contains(*lo_sup)) mk_bound_axiom(*a1, **lo_sup); + if (hi_inf1 != end && hi_inf != end && !visited.contains(*hi_inf)) mk_bound_axiom(*a1, **hi_inf); + if (hi_sup1 != end && hi_sup != end && !visited.contains(*hi_sup)) mk_bound_axiom(*a1, **hi_sup); + } + } + } + + struct compare_bounds { + bool operator()(lp::bound* a1, lp::bound* a2) const { return a1->get_value() < a2->get_value(); } + }; + + + lp_bounds::iterator first( + lp::bound_kind kind, + iterator it, + iterator end) { + for (; it != end; ++it) { + lp::bound* a = *it; + if (a->get_bound_kind() == kind) return it; + } + return end; + } + + lp_bounds::iterator next_inf( + lp::bound* a1, + lp::bound_kind kind, + iterator it, + iterator end, + bool& found_compatible) { + rational const & k1(a1->get_value()); + iterator result = end; + found_compatible = false; + for (; it != end; ++it) { + lp::bound * a2 = *it; + if (a1 == a2) continue; + if (a2->get_bound_kind() != kind) continue; + rational const & k2(a2->get_value()); + found_compatible = true; + if (k2 <= k1) { + result = it; + } + else { + break; + } + } + return result; + } + + lp_bounds::iterator next_sup( + lp::bound* a1, + lp::bound_kind kind, + iterator it, + iterator end, + bool& found_compatible) { + rational const & k1(a1->get_value()); + found_compatible = false; + for (; it != end; ++it) { + lp::bound * a2 = *it; + if (a1 == a2) continue; + if (a2->get_bound_kind() != kind) continue; + rational const & k2(a2->get_value()); + found_compatible = true; + if (k1 < k2) { + return it; + } + } + return end; + } + + void propagate_basic_bounds() { + for (auto const& bv : m_to_check) { + lp::bound& b = *m_bool_var2bound.find(bv); + propagate_bound(bv, ctx().get_assignment(bv) == l_true, b); + if (ctx().inconsistent()) break; + + } + m_to_check.reset(); + } + + // for glb lo': lo' < lo: + // lo <= x -> lo' <= x + // lo <= x -> ~(x <= lo') + // for lub hi': hi' > hi + // x <= hi -> x <= hi' + // x <= hi -> ~(x >= hi') + + void propagate_bound(bool_var bv, bool is_true, lp::bound& b) { + if (BP_NONE == propagation_mode()) { + return; + } + lp::bound_kind k = b.get_bound_kind(); + theory_var v = b.get_var(); + inf_rational val = b.get_value(is_true); + lp_bounds const& bounds = m_bounds[v]; + SASSERT(!bounds.empty()); + if (bounds.size() == 1) return; + if (m_unassigned_bounds[v] == 0) return; + + literal lit1(bv, !is_true); + literal lit2 = null_literal; + bool find_glb = (is_true == (k == lp::lower_t)); + if (find_glb) { + rational glb; + lp::bound* lb = 0; + for (unsigned i = 0; i < bounds.size(); ++i) { + lp::bound* b2 = bounds[i]; + if (b2 == &b) continue; + rational const& val2 = b2->get_value(); + if ((is_true ? val2 < val : val2 <= val) && (!lb || glb < val2)) { + lb = b2; + glb = val2; + } + } + if (!lb) return; + bool sign = lb->get_bound_kind() != lp::lower_t; + lit2 = literal(lb->get_bv(), sign); + } + else { + rational lub; + lp::bound* ub = 0; + for (unsigned i = 0; i < bounds.size(); ++i) { + lp::bound* b2 = bounds[i]; + if (b2 == &b) continue; + rational const& val2 = b2->get_value(); + if ((is_true ? val < val2 : val <= val2) && (!ub || val2 < lub)) { + ub = b2; + lub = val2; + } + } + if (!ub) return; + bool sign = ub->get_bound_kind() != lp::upper_t; + lit2 = literal(ub->get_bv(), sign); + } + TRACE("arith", + ctx().display_literal_verbose(tout, lit1); + ctx().display_literal_verbose(tout << " => ", lit2); + tout << "\n";); + updt_unassigned_bounds(v, -1); + ++m_stats.m_bound_propagations2; + m_params.reset(); + m_core.reset(); + m_eqs.reset(); + m_core.push_back(lit2); + m_params.push_back(parameter(symbol("farkas"))); + m_params.push_back(parameter(rational(1))); + m_params.push_back(parameter(rational(1))); + assign(lit2); + ++m_stats.m_bounds_propagations; + } + + void add_use_lists(lp::bound* b) { + theory_var v = b->get_var(); + lean::var_index vi = get_var_index(v); + if (m_solver->is_term(vi)) { + lean::lar_term const& term = m_solver->get_term(vi); + for (auto i = term.m_coeffs.begin(); i != term.m_coeffs.end(); ++i) { + lean::var_index wi = i->first; + unsigned w = m_var_index2theory_var[wi]; + m_use_list.reserve(w + 1, ptr_vector()); + m_use_list[w].push_back(b); + } + } + } + + void del_use_lists(lp::bound* b) { + theory_var v = b->get_var(); + lean::var_index vi = m_theory_var2var_index[v]; + if (m_solver->is_term(vi)) { + lean::lar_term const& term = m_solver->get_term(vi); + for (auto i = term.m_coeffs.begin(); i != term.m_coeffs.end(); ++i) { + lean::var_index wi = i->first; + unsigned w = m_var_index2theory_var[wi]; + SASSERT(m_use_list[w].back() == b); + m_use_list[w].pop_back(); + } + } + } + + // + // propagate bounds to compound terms + // The idea is that if bounds on all variables in an inequality ax + by + cz >= k + // have been assigned we may know the truth value of the inequality by using simple + // bounds propagation. + // + void propagate_bound_compound(bool_var bv, bool is_true, lp::bound& b) { + lp::bound_kind k = b.get_bound_kind(); + theory_var v = b.get_var(); + TRACE("arith", tout << mk_pp(get_owner(v), m) << "\n";); + if (static_cast(v) >= m_use_list.size()) { + return; + } + for (auto const& vb : m_use_list[v]) { + if (ctx().get_assignment(vb->get_bv()) != l_undef) { + TRACE("arith_verbose", display_bound(tout << "assigned ", *vb) << "\n";); + continue; + } + inf_rational r; + // x + y + // x >= 0, y >= 1 -> x + y >= 1 + // x <= 0, y <= 2 -> x + y <= 2 + literal lit = null_literal; + if (lp::lower_t == vb->get_bound_kind()) { + if (get_glb(*vb, r) && r >= vb->get_value()) { // vb is assigned true + lit = literal(vb->get_bv(), false); + } + else if (get_lub(*vb, r) && r < vb->get_value()) { // vb is assigned false + lit = literal(vb->get_bv(), true); + } + } + else { + if (get_glb(*vb, r) && r > vb->get_value()) { // VB <= value < val(VB) + lit = literal(vb->get_bv(), true); + } + else if (get_lub(*vb, r) && r <= vb->get_value()) { // val(VB) <= value + lit = literal(vb->get_bv(), false); + } + } + + if (lit != null_literal) { + TRACE("arith", + ctx().display_literals_verbose(tout, m_core); + tout << "\n --> "; + ctx().display_literal_verbose(tout, lit); + tout << "\n"; + ); + + + assign(lit); + } + else { + TRACE("arith_verbose", display_bound(tout << "skip ", *vb) << "\n";); + } + } + } + + bool get_lub(lp::bound const& b, inf_rational& lub) { + return get_bound(b, lub, true); + } + + bool get_glb(lp::bound const& b, inf_rational& glb) { + return get_bound(b, glb, false); + } + + std::ostream& display_bound(std::ostream& out, lp::bound const& b) { + return out << mk_pp(ctx().bool_var2expr(b.get_bv()), m); + } + + bool get_bound(lp::bound const& b, inf_rational& r, bool is_lub) { + m_core.reset(); + m_eqs.reset(); + m_params.reset(); + r.reset(); + theory_var v = b.get_var(); + lean::var_index vi = m_theory_var2var_index[v]; + SASSERT(m_solver->is_term(vi)); + lean::lar_term const& term = m_solver->get_term(vi); + for (auto const coeff : term.m_coeffs) { + lean::var_index wi = coeff.first; + lean::constraint_index ci; + rational value; + bool is_strict; + if (coeff.second.is_neg() == is_lub) { + // -3*x ... <= lub based on lower bound for x. + if (!m_solver->has_lower_bound(wi, ci, value, is_strict)) { + return false; + } + if (is_strict) { + r += inf_rational(rational::zero(), coeff.second.is_pos()); + } + } + else { + if (!m_solver->has_upper_bound(wi, ci, value, is_strict)) { + return false; + } + if (is_strict) { + r += inf_rational(rational::zero(), coeff.second.is_pos()); + } + } + r += value * coeff.second; + set_evidence(ci); + } + TRACE("arith_verbose", tout << (is_lub?"lub":"glb") << " is " << r << "\n";); + return true; + } + + void assert_bound(bool_var bv, bool is_true, lp::bound& b) { + if (m_solver->get_status() == lean::lp_status::INFEASIBLE) { + return; + } + scoped_internalize_state st(*this); + st.vars().push_back(b.get_var()); + st.coeffs().push_back(rational::one()); + init_left_side(st); + lean::lconstraint_kind k = lean::EQ; + switch (b.get_bound_kind()) { + case lp::lower_t: + k = is_true ? lean::GE : lean::LT; + break; + case lp::upper_t: + k = is_true ? lean::LE : lean::GT; + break; + } + if (k == lean::LT || k == lean::LE) { + ++m_stats.m_assert_lower; + } + else { + ++m_stats.m_assert_upper; + } + auto vi = get_var_index(b.get_var()); + auto ci = m_solver->add_var_bound(vi, k, b.get_value()); + TRACE("arith", tout << "v" << b.get_var() << "\n";); + add_ineq_constraint(ci, literal(bv, !is_true)); + + propagate_eqs(vi, ci, k, b); + } + + // + // fixed equalities. + // A fixed equality is inferred if there are two variables v1, v2 whose + // upper and lower bounds coincide. + // Then the equality v1 == v2 is propagated to the core. + // + + typedef std::pair constraint_bound; + vector m_lower_terms; + vector m_upper_terms; + typedef std::pair value_sort_pair; + typedef pair_hash, bool_hash> value_sort_pair_hash; + typedef map > value2var; + value2var m_fixed_var_table; + const lean::constraint_index null_index = static_cast(-1); + + void propagate_eqs(lean::var_index vi, lean::constraint_index ci, lean::lconstraint_kind k, lp::bound& b) { + if (propagate_eqs()) { + rational const& value = b.get_value(); + if (k == lean::GE) { + set_lower_bound(vi, ci, value); + if (has_upper_bound(vi, ci, value)) { + fixed_var_eh(b.get_var(), value); + } + } + else if (k == lean::LE) { + set_upper_bound(vi, ci, value); + if (has_lower_bound(vi, ci, value)) { + fixed_var_eh(b.get_var(), value); + } + } + } + } + + bool dump_lemmas() const { return m_arith_params.m_arith_dump_lemmas; } + + bool propagate_eqs() const { return m_arith_params.m_arith_propagate_eqs && m_num_conflicts < m_arith_params.m_arith_propagation_threshold; } + + bound_prop_mode propagation_mode() const { return m_num_conflicts < m_arith_params.m_arith_propagation_threshold ? m_arith_params.m_arith_bound_prop : BP_NONE; } + + unsigned small_lemma_size() const { return m_arith_params.m_arith_small_lemma_size; } + + bool proofs_enabled() const { return m.proofs_enabled(); } + + bool use_tableau() const { return m_lp_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); } + + void set_lower_bound(lean::var_index vi, lean::constraint_index ci, rational const& v) { set_bound(vi, ci, v, true); } + + void set_bound(lean::var_index vi, lean::constraint_index ci, rational const& v, bool is_lower) { + if (!m_solver->is_term(vi)) { + // m_solver already tracks bounds on proper variables, but not on terms. + return; + } + 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())); + } + constraint_bound& b = vec[ti]; + if (b.first == null_index || (is_lower? b.second < v : b.second > v)) { + ctx().push_trail(vector_value_trail(vec, ti)); + b.first = ci; + b.second = v; + } + } + + bool has_upper_bound(lean::var_index vi, lean::constraint_index& ci, rational const& bound) { return has_bound(vi, ci, bound, false); } + + bool has_lower_bound(lean::var_index vi, lean::constraint_index& ci, rational const& bound) { return has_bound(vi, ci, bound, true); } + + bool has_bound(lean::var_index vi, lean::constraint_index& ci, rational const& bound, bool is_lower) { + + if (m_solver->is_term(vi)) { + + lean::var_index ti = m_solver->adjust_term_index(vi); + theory_var v = m_term_index2theory_var.get(ti, null_theory_var); + 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; + return bound == val; + } + + auto& vec = is_lower ? m_lower_terms : m_upper_terms; + if (vec.size() > ti) { + constraint_bound& b = vec[ti]; + ci = b.first; + return ci != null_index && bound == b.second; + } + else { + return false; + + } + } + else { + bool is_strict = false; + rational b; + if (is_lower) { + return m_solver->has_lower_bound(vi, ci, b, is_strict) && b == bound && !is_strict; + } + else { + return m_solver->has_upper_bound(vi, ci, b, is_strict) && b == bound && !is_strict; + } + } + } + + + bool is_equal(theory_var x, theory_var y) const { return get_enode(x)->get_root() == get_enode(y)->get_root(); } + + + void fixed_var_eh(theory_var v1, rational const& bound) { + theory_var v2; + value_sort_pair key(bound, is_int(v1)); + if (m_fixed_var_table.find(key, v2)) { + if (static_cast(v2) < th.get_num_vars() && !is_equal(v1, v2)) { + auto vi1 = get_var_index(v1); + auto vi2 = get_var_index(v2); + lean::constraint_index ci1, ci2, ci3, ci4; + TRACE("arith", tout << "fixed: " << mk_pp(get_owner(v1), m) << " " << mk_pp(get_owner(v2), m) << " " << bound << " " << has_lower_bound(vi2, ci3, bound) << "\n";); + if (has_lower_bound(vi2, ci3, bound) && has_upper_bound(vi2, ci4, bound)) { + VERIFY (has_lower_bound(vi1, ci1, bound)); + VERIFY (has_upper_bound(vi1, ci2, bound)); + ++m_stats.m_fixed_eqs; + m_core.reset(); + m_eqs.reset(); + set_evidence(ci1); + set_evidence(ci2); + set_evidence(ci3); + set_evidence(ci4); + enode* x = get_enode(v1); + enode* y = get_enode(v2); + justification* js = + ctx().mk_justification( + ext_theory_eq_propagation_justification( + get_id(), ctx().get_region(), m_core.size(), m_core.c_ptr(), m_eqs.size(), m_eqs.c_ptr(), x, y, 0, 0)); + + TRACE("arith", + for (unsigned i = 0; i < m_core.size(); ++i) { + ctx().display_detailed_literal(tout, m_core[i]); + tout << "\n"; + } + for (unsigned i = 0; i < m_eqs.size(); ++i) { + tout << mk_pp(m_eqs[i].first->get_owner(), m) << " = " << mk_pp(m_eqs[i].second->get_owner(), m) << "\n"; + } + tout << " ==> "; + tout << mk_pp(x->get_owner(), m) << " = " << mk_pp(y->get_owner(), m) << "\n"; + ); + + // parameters are TBD. + SASSERT(validate_eq(x, y)); + ctx().assign_eq(x, y, eq_justification(js)); + } + } + else { + // bounds on v2 were changed. + m_fixed_var_table.insert(key, v1); + } + } + else { + m_fixed_var_table.insert(key, v1); + } + } + + lbool make_feasible() { + reset_variable_values(); + ++m_stats.m_make_feasible; + if (m_solver->A_r().column_count() > m_stats.m_max_cols) + m_stats.m_max_cols = m_solver->A_r().column_count(); + 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: + return l_false; + case lean::lp_status::FEASIBLE: + case lean::lp_status::OPTIMAL: + // SASSERT(m_solver->all_constraints_hold()); + return l_true; + case lean::lp_status::TIME_EXHAUSTED: + + default: + TRACE("arith", tout << "status treated as inconclusive: " << status << "\n";); + // TENTATIVE_UNBOUNDED, UNBOUNDED, TENTATIVE_DUAL_UNBOUNDED, DUAL_UNBOUNDED, + // FLOATING_POINT_ERROR, TIME_EXAUSTED, ITERATIONS_EXHAUSTED, EMPTY, UNSTABLE + return l_undef; + } + } + + vector> m_explanation; + literal_vector m_core; + svector m_eqs; + vector m_params; + lean::constraint_index const null_constraint_index = UINT_MAX; + + void set_evidence(lean::constraint_index idx) { + if (idx == null_constraint_index) { + return; + } + switch (m_constraint_sources[idx]) { + case inequality_source: { + literal lit = m_inequalities[idx]; + SASSERT(lit != null_literal); + m_core.push_back(lit); + break; + } + case equality_source: { + SASSERT(m_equalities[idx].first != nullptr); + SASSERT(m_equalities[idx].second != nullptr); + m_eqs.push_back(m_equalities[idx]); + break; + } + case definition_source: { + // skip definitions (these are treated as hard constraints) + break; + } + default: + UNREACHABLE(); + break; + } + } + + void set_conflict() { + m_eqs.reset(); + m_core.reset(); + m_params.reset(); + m_explanation.clear(); + m_solver->get_infeasibility_explanation(m_explanation); + // m_solver->shrink_explanation_to_minimum(m_explanation); // todo, enable when perf is fixed + /* + static unsigned cn = 0; + static unsigned num_l = 0; + num_l+=m_explanation.size(); + std::cout << num_l / (++cn) << "\n"; + */ + ++m_num_conflicts; + ++m_stats.m_conflicts; + TRACE("arith", tout << "scope: " << ctx().get_scope_level() << "\n"; display_evidence(tout, m_explanation); ); + TRACE("arith", display(tout);); + for (auto const& ev : m_explanation) { + if (!ev.first.is_zero()) { + set_evidence(ev.second); + } + } + SASSERT(validate_conflict()); + ctx().set_conflict( + ctx().mk_justification( + ext_theory_conflict_justification( + get_id(), ctx().get_region(), + m_core.size(), m_core.c_ptr(), + m_eqs.size(), m_eqs.c_ptr(), m_params.size(), m_params.c_ptr()))); + } + + justification * why_is_diseq(theory_var v1, theory_var v2) { + return 0; + } + + void reset_eh() { + m_arith_eq_adapter.reset_eh(); + m_solver = 0; + m_not_handled = nullptr; + del_bounds(0); + m_unassigned_bounds.reset(); + m_asserted_qhead = 0; + m_scopes.reset(); + m_stats.reset(); + m_to_check.reset(); + } + + void init_model(model_generator & mg) { + init_variable_values(); + m_factory = alloc(arith_factory, m); + mg.register_factory(m_factory); + TRACE("arith", display(tout);); + } + + model_value_proc * mk_value(enode * n, model_generator & mg) { + theory_var v = n->get_th_var(get_id()); + return alloc(expr_wrapper_proc, m_factory->mk_value(get_value(v), m.get_sort(n->get_owner()))); + } + + bool get_value(enode* n, expr_ref& r) { + theory_var v = n->get_th_var(get_id()); + if (can_get_value(v)) { + r = a.mk_numeral(get_value(v), is_int(n)); + return true; + } + else { + return false; + } + } + + bool validate_eq_in_model(theory_var v1, theory_var v2, bool is_true) const { + SASSERT(v1 != null_theory_var); + SASSERT(v2 != null_theory_var); + return (get_value(v1) == get_value(v2)) == is_true; + } + + // Auxiliary verification utilities. + + bool validate_conflict() { + if (dump_lemmas()) { + ctx().display_lemma_as_smt_problem(m_core.size(), m_core.c_ptr(), m_eqs.size(), m_eqs.c_ptr(), false_literal); + } + context nctx(m, ctx().get_fparams(), ctx().get_params()); + add_background(nctx); + bool result = l_true != nctx.check(); + CTRACE("arith", !result, ctx().display_lemma_as_smt_problem(tout, m_core.size(), m_core.c_ptr(), m_eqs.size(), m_eqs.c_ptr(), false_literal); + display(tout);); + return result; + } + + bool validate_assign(literal lit) { + if (dump_lemmas()) { + ctx().display_lemma_as_smt_problem(m_core.size(), m_core.c_ptr(), m_eqs.size(), m_eqs.c_ptr(), lit); + } + context nctx(m, ctx().get_fparams(), ctx().get_params()); + m_core.push_back(~lit); + add_background(nctx); + m_core.pop_back(); + bool result = l_true != nctx.check(); + CTRACE("arith", !result, ctx().display_lemma_as_smt_problem(tout, m_core.size(), m_core.c_ptr(), m_eqs.size(), m_eqs.c_ptr(), lit); + display(tout);); + return result; + } + + bool validate_eq(enode* x, enode* y) { + context nctx(m, ctx().get_fparams(), ctx().get_params()); + add_background(nctx); + nctx.assert_expr(m.mk_not(m.mk_eq(x->get_owner(), y->get_owner()))); + return l_true != nctx.check(); + } + + void add_background(context& nctx) { + for (unsigned i = 0; i < m_core.size(); ++i) { + expr_ref tmp(m); + ctx().literal2expr(m_core[i], tmp); + nctx.assert_expr(tmp); + } + for (unsigned i = 0; i < m_eqs.size(); ++i) { + nctx.assert_expr(m.mk_eq(m_eqs[i].first->get_owner(), m_eqs[i].second->get_owner())); + } + } + + theory_lra::inf_eps value(theory_var v) { + lean::impq ival = get_ivalue(v); + return inf_eps(0, inf_rational(ival.x, ival.y)); + } + + theory_lra::inf_eps maximize(theory_var v, expr_ref& blocker, bool& has_shared) { + lean::var_index vi = m_theory_var2var_index.get(v, UINT_MAX); + vector > coeffs; + rational coeff; + if (m_solver->is_term(vi)) { + const lean::lar_term& term = m_solver->get_term(vi); + for (auto & ti : term.m_coeffs) { + coeffs.push_back(std::make_pair(ti.second, ti.first)); + } + coeff = term.m_v; + } + else { + coeffs.push_back(std::make_pair(rational::one(), vi)); + coeff = rational::zero(); + } + lean::impq term_max; + if (m_solver->maximize_term(coeffs, term_max)) { + blocker = mk_gt(v); + inf_rational val(term_max.x + coeff, term_max.y); + return inf_eps(rational::zero(), val); + } + else { + TRACE("arith", tout << "Unbounded " << v << "\n";); + has_shared = false; + blocker = m.mk_false(); + return inf_eps(rational::one(), inf_rational()); + } + } + + expr_ref mk_gt(theory_var v) { + lean::impq val = get_ivalue(v); + expr* obj = get_enode(v)->get_owner(); + rational r = val.x; + expr_ref e(m); + if (a.is_int(m.get_sort(obj))) { + if (r.is_int()) { + r += rational::one(); + } + else { + r = ceil(r); + } + e = a.mk_numeral(r, m.get_sort(obj)); + e = a.mk_ge(obj, e); + } + else { + e = a.mk_numeral(r, m.get_sort(obj)); + if (val.y.is_neg()) { + e = a.mk_ge(obj, e); + } + else { + e = a.mk_gt(obj, e); + } + } + TRACE("opt", tout << "v" << v << " " << val << " " << r << " " << e << "\n";); + return e; + } + + theory_var add_objective(app* term) { + return internalize_def(term); + } + + app_ref mk_obj(theory_var v) { + lean::var_index vi = m_theory_var2var_index[v]; + if (m_solver->is_term(vi)) { + expr_ref_vector args(m); + const lean::lar_term& term = m_solver->get_term(vi); + for (auto & ti : term.m_coeffs) { + theory_var w = m_var_index2theory_var[ti.first]; + expr* o = get_enode(w)->get_owner(); + args.push_back(a.mk_mul(a.mk_numeral(ti.second, a.is_int(o)), o)); + } + rational r = term.m_v; + args.push_back(a.mk_numeral(r, r.is_int())); + return app_ref(a.mk_add(args.size(), args.c_ptr()), m); + } + else { + theory_var w = m_var_index2theory_var[vi]; + return app_ref(get_enode(w)->get_owner(), m); + } + } + + expr_ref mk_ge(filter_model_converter& fm, theory_var v, inf_rational const& val) { + rational r = val.get_rational(); + bool is_strict = val.get_infinitesimal().is_pos(); + app_ref b(m); + if (is_strict) { + b = a.mk_le(mk_obj(v), a.mk_numeral(r, r.is_int())); + } + else { + b = a.mk_ge(mk_obj(v), a.mk_numeral(r, r.is_int())); + } + if (!ctx().b_internalized(b)) { + fm.insert(b->get_decl()); + bool_var bv = ctx().mk_bool_var(b); + ctx().set_var_theory(bv, get_id()); + // ctx().set_enode_flag(bv, true); + lp::bound_kind bkind = lp::bound_kind::lower_t; + if (is_strict) bkind = lp::bound_kind::upper_t; + lp::bound* a = alloc(lp::bound, bv, v, r, bkind); + mk_bound_axioms(*a); + updt_unassigned_bounds(v, +1); + m_bounds[v].push_back(a); + m_bounds_trail.push_back(v); + m_bool_var2bound.insert(bv, a); + TRACE("arith", tout << mk_pp(b, m) << "\n";); + } + if (is_strict) { + b = m.mk_not(b); + } + TRACE("arith", tout << b << "\n";); + return expr_ref(b, m); + + } + + + void display(std::ostream & out) const { + if (m_solver) { + m_solver->print_constraints(out); + m_solver->print_terms(out); + } + unsigned nv = th.get_num_vars(); + for (unsigned v = 0; v < nv; ++v) { + out << "v" << v; + if (can_get_value(v)) out << ", value: " << get_value(v); + out << ", shared: " << ctx().is_shared(get_enode(v)) + << ", rel: " << ctx().is_relevant(get_enode(v)) + << ", def: "; th.display_var_flat_def(out, v) << "\n"; + } + } + + void display_evidence(std::ostream& out, vector> const& evidence) { + for (auto const& ev : evidence) { + expr_ref e(m); + SASSERT(!ev.first.is_zero()); + if (ev.first.is_zero()) { + continue; + } + unsigned idx = ev.second; + switch (m_constraint_sources.get(idx, null_source)) { + case inequality_source: { + literal lit = m_inequalities[idx]; + ctx().literal2expr(lit, e); + out << e << " " << ctx().get_assignment(lit) << "\n"; + break; + } + case equality_source: + out << mk_pp(m_equalities[idx].first->get_owner(), m) << " = " + << mk_pp(m_equalities[idx].second->get_owner(), m) << "\n"; + break; + case definition_source: { + theory_var v = m_definitions[idx]; + out << "def: v" << v << " := " << mk_pp(th.get_enode(v)->get_owner(), m) << "\n"; + break; + } + case null_source: + default: + UNREACHABLE(); + break; + } + } + for (auto const& ev : evidence) { + m_solver->print_constraint(ev.second, out << ev.first << ": "); + } + } + + void collect_statistics(::statistics & st) const { + m_arith_eq_adapter.collect_statistics(st); + st.update("arith-lower", m_stats.m_assert_lower); + st.update("arith-upper", m_stats.m_assert_upper); + st.update("arith-rows", m_stats.m_add_rows); + st.update("arith-propagations", m_stats.m_bounds_propagations); + st.update("arith-iterations", m_stats.m_num_iterations); + st.update("arith-factorizations", m_stats.m_num_factorizations); + st.update("arith-pivots", m_stats.m_need_to_solve_inf); + st.update("arith-plateau-iterations", m_stats.m_num_iterations_with_no_progress); + st.update("arith-fixed-eqs", m_stats.m_fixed_eqs); + st.update("arith-conflicts", m_stats.m_conflicts); + st.update("arith-bound-propagations-lp", m_stats.m_bound_propagations1); + st.update("arith-bound-propagations-cheap", m_stats.m_bound_propagations2); + st.update("arith-diseq", m_stats.m_assert_diseq); + st.update("arith-make-feasible", m_stats.m_make_feasible); + st.update("arith-max-columns", m_stats.m_max_cols); + st.update("arith-max-rows", m_stats.m_max_rows); + } + }; + + theory_lra::theory_lra(ast_manager& m, theory_arith_params& p): + theory(m.get_family_id("arith")) { + m_imp = alloc(imp, *this, m, p); + } + theory_lra::~theory_lra() { + dealloc(m_imp); + } + theory* theory_lra::mk_fresh(context* new_ctx) { + return alloc(theory_lra, new_ctx->get_manager(), new_ctx->get_fparams()); + } + void theory_lra::init(context * ctx) { + theory::init(ctx); + } + bool theory_lra::internalize_atom(app * atom, bool gate_ctx) { + return m_imp->internalize_atom(atom, gate_ctx); + } + bool theory_lra::internalize_term(app * term) { + return m_imp->internalize_term(term); + } + void theory_lra::internalize_eq_eh(app * atom, bool_var v) { + m_imp->internalize_eq_eh(atom, v); + } + void theory_lra::assign_eh(bool_var v, bool is_true) { + m_imp->assign_eh(v, is_true); + } + void theory_lra::new_eq_eh(theory_var v1, theory_var v2) { + m_imp->new_eq_eh(v1, v2); + } + bool theory_lra::use_diseqs() const { + return m_imp->use_diseqs(); + } + void theory_lra::new_diseq_eh(theory_var v1, theory_var v2) { + m_imp->new_diseq_eh(v1, v2); + } + void theory_lra::push_scope_eh() { + theory::push_scope_eh(); + m_imp->push_scope_eh(); + } + void theory_lra::pop_scope_eh(unsigned num_scopes) { + m_imp->pop_scope_eh(num_scopes); + theory::pop_scope_eh(num_scopes); + } + void theory_lra::restart_eh() { + m_imp->restart_eh(); + } + void theory_lra::relevant_eh(app* e) { + m_imp->relevant_eh(e); + } + void theory_lra::init_search_eh() { + m_imp->init_search_eh(); + } + final_check_status theory_lra::final_check_eh() { + return m_imp->final_check_eh(); + } + bool theory_lra::is_shared(theory_var v) const { + return m_imp->is_shared(v); + } + bool theory_lra::can_propagate() { + return m_imp->can_propagate(); + } + void theory_lra::propagate() { + m_imp->propagate(); + } + justification * theory_lra::why_is_diseq(theory_var v1, theory_var v2) { + return m_imp->why_is_diseq(v1, v2); + } + void theory_lra::reset_eh() { + m_imp->reset_eh(); + } + void theory_lra::init_model(model_generator & m) { + m_imp->init_model(m); + } + model_value_proc * theory_lra::mk_value(enode * n, model_generator & mg) { + return m_imp->mk_value(n, mg); + } + bool theory_lra::get_value(enode* n, expr_ref& r) { + return m_imp->get_value(n, r); + } + bool theory_lra::validate_eq_in_model(theory_var v1, theory_var v2, bool is_true) const { + return m_imp->validate_eq_in_model(v1, v2, is_true); + } + void theory_lra::display(std::ostream & out) const { + m_imp->display(out); + } + void theory_lra::collect_statistics(::statistics & st) const { + m_imp->collect_statistics(st); + } + theory_lra::inf_eps theory_lra::value(theory_var v) { + return m_imp->value(v); + } + theory_lra::inf_eps theory_lra::maximize(theory_var v, expr_ref& blocker, bool& has_shared) { + return m_imp->maximize(v, blocker, has_shared); + } + theory_var theory_lra::add_objective(app* term) { + return m_imp->add_objective(term); + } + expr_ref theory_lra::mk_ge(filter_model_converter& fm, theory_var v, inf_rational const& val) { + return m_imp->mk_ge(fm, v, val); + } + + + +} diff --git a/src/smt/theory_lra.h b/src/smt/theory_lra.h new file mode 100644 index 000000000..c91363848 --- /dev/null +++ b/src/smt/theory_lra.h @@ -0,0 +1,96 @@ +/*++ +Copyright (c) 2016 Microsoft Corporation + +Module Name: + + theory_lra.h + +Abstract: + + + +Author: + + Lev Nachmanson (levnach) 2016-25-3 + Nikolaj Bjorner (nbjorner) + +Revision History: + + +--*/ +#pragma once + +#include "theory_opt.h" + +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); + virtual ~theory_lra(); + virtual theory* mk_fresh(context* new_ctx); + virtual char const* get_name() const { return "lra"; } + + virtual void init(context * ctx); + + virtual bool internalize_atom(app * atom, bool gate_ctx); + + virtual bool internalize_term(app * term); + + virtual void internalize_eq_eh(app * atom, bool_var v); + + virtual void assign_eh(bool_var v, bool is_true); + + virtual void new_eq_eh(theory_var v1, theory_var v2); + + virtual bool use_diseqs() const; + + virtual void new_diseq_eh(theory_var v1, theory_var v2); + + virtual void push_scope_eh(); + + virtual void pop_scope_eh(unsigned num_scopes); + + virtual void restart_eh(); + + virtual void relevant_eh(app* e); + + virtual void init_search_eh(); + + virtual final_check_status final_check_eh(); + + virtual bool is_shared(theory_var v) const; + + virtual bool can_propagate(); + + virtual void propagate(); + + virtual justification * why_is_diseq(theory_var v1, theory_var v2); + + // virtual void flush_eh(); + + virtual void reset_eh(); + + virtual void init_model(model_generator & m); + + virtual model_value_proc * mk_value(enode * n, model_generator & mg); + + virtual bool get_value(enode* n, expr_ref& r); + + virtual bool validate_eq_in_model(theory_var v1, theory_var v2, bool is_true) const; + + virtual void display(std::ostream & out) const; + + virtual void collect_statistics(::statistics & st) const; + + // optimization + virtual inf_eps value(theory_var); + virtual inf_eps maximize(theory_var v, expr_ref& blocker, bool& has_shared); + virtual theory_var add_objective(app* term); + virtual expr_ref mk_ge(filter_model_converter& fm, theory_var v, inf_rational const& val); + + + }; + +} diff --git a/src/test/argument_parser.h b/src/test/argument_parser.h new file mode 100644 index 000000000..706167f49 --- /dev/null +++ b/src/test/argument_parser.h @@ -0,0 +1,144 @@ +/* +Copyright (c) 2013 Microsoft Corporation. All rights reserved. +Released under Apache 2.0 license as described in the file LICENSE. + +Author: Lev Nachmanson +*/ + +#include +#include +#include +#include +#include + +namespace lean { +class argument_parser { + std::unordered_map m_options; + std::unordered_map m_options_with_after_string; + std::set m_used_options; + std::unordered_map m_used_options_with_after_string; + std::vector m_free_args; + std::vector m_args; + +public: + std::string m_error_message; + argument_parser(unsigned argn, char * const* args) { + for (unsigned i = 0; i < argn; i++) { + m_args.push_back(std::string(args[i])); + } + } + + void add_option(std::string s) { + add_option_with_help_string(s, ""); + } + + void add_option_with_help_string(std::string s, std::string help_string) { + m_options[s]=help_string; + } + + void add_option_with_after_string(std::string s) { + add_option_with_after_string_with_help(s, ""); + } + + void add_option_with_after_string_with_help(std::string s, std::string help_string) { + m_options_with_after_string[s]=help_string; + } + + + bool parse() { + bool status_is_ok = true; + for (unsigned i = 0; i < m_args.size(); i++) { + std::string ar = m_args[i]; + if (m_options.find(ar) != m_options.end() ) + m_used_options.insert(ar); + else if (m_options_with_after_string.find(ar) != m_options_with_after_string.end()) { + if (i == m_args.size() - 1) { + m_error_message = "Argument is missing after "+ar; + return false; + } + i++; + m_used_options_with_after_string[ar] = m_args[i]; + } else { + if (starts_with(ar, "-") || starts_with(ar, "//")) + status_is_ok = false; + + m_free_args.push_back(ar); + } + } + return status_is_ok; + } + + bool contains(std::unordered_map & m, std::string s) { + return m.find(s) != m.end(); + } + + bool contains(std::set & m, std::string s) { + return m.find(s) != m.end(); + } + + bool option_is_used(std::string option) { + return contains(m_used_options, option) || contains(m_used_options_with_after_string, option); + } + + std::string get_option_value(std::string option) { + auto t = m_used_options_with_after_string.find(option); + if (t != m_used_options_with_after_string.end()){ + return t->second; + } + return std::string(); + } + + bool starts_with(std::string s, char const * prefix) { + return starts_with(s, std::string(prefix)); + } + + bool starts_with(std::string s, std::string prefix) { + return s.substr(0, prefix.size()) == prefix; + } + + std::string usage_string() { + std::string ret = ""; + std::vector unknown_options; + for (auto t : m_free_args) { + if (starts_with(t, "-") || starts_with(t, "\\")) { + unknown_options.push_back(t); + } + } + if (unknown_options.size()) { + ret = "Unknown options:"; + } + for (auto unknownOption : unknown_options) { + ret += unknownOption; + ret += ","; + } + ret += "\n"; + ret += "Usage:\n"; + for (auto allowed_option : m_options) + ret += allowed_option.first + " " + (allowed_option.second.size() == 0 ? std::string("") : std::string("/") + allowed_option.second) + std::string("\n"); + for (auto s : m_options_with_after_string) { + ret += s.first + " " + (s.second.size() == 0? " \"option value\"":("\""+ s.second+"\"")) + "\n"; + } + return ret; + } + + void print() { + if (m_used_options.size() == 0 && m_used_options_with_after_string.size() == 0 && m_free_args.size() == 0) { + std::cout << "no options are given" << std::endl; + return; + } + std::cout << "options are: " << std::endl; + for (std::string s : m_used_options) { + std::cout << s << std::endl; + } + for (auto & t : m_used_options_with_after_string) { + std::cout << t.first << " " << t.second << std::endl; + } + if (m_free_args.size() > 0) { + std::cout << "free arguments are: " << std::endl; + for (auto & t : m_free_args) { + std::cout << t << " " << std::endl; + } + } + } +}; +} diff --git a/src/test/lp.cpp b/src/test/lp.cpp new file mode 100644 index 000000000..7829d81f5 --- /dev/null +++ b/src/test/lp.cpp @@ -0,0 +1,3232 @@ +/* +Copyright (c) 2017 Microsoft Corporation. All rights reserved. +Author: Lev Nachmanson +*/ +#include +#if _LINUX_ +#include +#endif +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include "util/lp/lp_utils.h" +#include "util/lp/lp_primal_simplex.h" +#include "util/lp/mps_reader.h" +#include "smt_reader.h" +#include "util/lp/binary_heap_priority_queue.h" +#include "argument_parser.h" +#include "test_file_reader.h" +#include "util/lp/indexed_value.h" +#include "util/lp/lar_solver.h" +#include "util/lp/numeric_pair.h" +#include "util/lp/binary_heap_upair_queue.h" +#include "util/lp/stacked_value.h" +#include "util/lp/stacked_unordered_set.h" +#include "util/lp/int_set.h" +#include "util/stopwatch.h" +namespace lean { +unsigned seed = 1; + +struct simple_column_namer:public column_namer +{ + std::string get_column_name(unsigned j) const override { + return std::string("x") + T_to_string(j); + } +}; + + +template +void test_matrix(sparse_matrix & a) { + auto m = a.dimension(); + +// copy a to b in the reversed order + sparse_matrix b(m); + std::cout << "copy b to a"<< std::endl; + for (int row = m - 1; row >= 0; row--) + for (int col = m - 1; col >= 0; col --) { + b(row, col) = (T const&) a(row, col); + } + + + std::cout << "zeroing b in the reverse order"<< std::endl; + for (int row = m - 1; row >= 0; row--) + for (int col = m - 1; col >= 0; col --) + b.set(row, col, T(0)); + + + + for (unsigned row = 0; row < m; row ++) + for (unsigned col = 0; col < m; col ++) + a.set(row, col, T(0)); + + + unsigned i = my_random() % m; + unsigned j = my_random() % m; + + auto t = T(1); + + a.set(i, j, t); + + lean_assert(a.get(i, j) == t); + + unsigned j1; + if (j < m - 1) { + j1 = m - 1; + a.set(i, j1, T(2)); + } +} + +void tst1() { + std::cout << "testing the minimial matrix with 1 row and 1 column" << std::endl; + sparse_matrix m0(1); + m0.set(0, 0, 1); + // print_matrix(m0); + m0.set(0, 0, 0); + // print_matrix(m0); + test_matrix(m0); + + unsigned rows = 2; + sparse_matrix m(rows); + std::cout << "setting m(0,1)=" << std::endl; + + m.set(0, 1, 11); + m.set(0, 0, 12); + + // print_matrix(m); + + test_matrix(m); + + sparse_matrix m1(2); + m1.set(0, 0, 2); + m1.set(1, 0, 3); + // print_matrix(m1); + std::cout << " zeroing matrix 2 by 2" << std::endl; + m1.set(0, 0, 0); + m1.set(1, 0, 0); + // print_matrix(m1); + + test_matrix(m1); + + + std::cout << "printing zero matrix 3 by 1" << std::endl; + sparse_matrix m2(3); + // print_matrix(m2); + + m2.set(0, 0, 1); + m2.set(2, 0, 2); + std::cout << "printing matrix 3 by 1 with a gap" << std::endl; + // print_matrix(m2); + + test_matrix(m2); + + sparse_matrix m10by9(10); + m10by9.set(0, 1, 1); + + m10by9(0, 1) = 4; + + double test = m10by9(0, 1); + + std::cout << "got " << test << std::endl; + + + m10by9.set(0, 8, 8); + m10by9.set(3, 4, 7); + m10by9.set(3, 2, 5); + m10by9.set(3, 8, 99); + m10by9.set(3, 2, 6); + m10by9.set(1, 8, 9); + m10by9.set(4, 0, 40); + m10by9.set(0, 0, 10); + + std::cout << "printing matrix 10 by 9" << std::endl; + // print_matrix(m10by9); + + + test_matrix(m10by9); + std::cout <<"zeroing m10by9\n"; +#ifdef LEAN_DEBUG + for (unsigned int i = 0; i < m10by9.dimension(); i++) + for (unsigned int j = 0; j < m10by9.column_count(); j++) + m10by9.set(i, j, 0); +#endif + // print_matrix(m10by9); +} + +vector allocate_basis_heading(unsigned count) { // the rest of initilization will be handled by lu_QR + vector basis_heading(count, -1); + return basis_heading; +} + + +void init_basic_part_of_basis_heading(vector & basis, vector & basis_heading) { + lean_assert(basis_heading.size() >= basis.size()); + unsigned m = basis.size(); + for (unsigned i = 0; i < m; i++) { + unsigned column = basis[i]; + basis_heading[column] = i; + } +} + +void init_non_basic_part_of_basis_heading(vector & basis_heading, vector & non_basic_columns) { + non_basic_columns.clear(); + for (int j = basis_heading.size(); j--;){ + if (basis_heading[j] < 0) { + non_basic_columns.push_back(j); + // the index of column j in m_nbasis is (- basis_heading[j] - 1) + basis_heading[j] = - static_cast(non_basic_columns.size()); + } + } +} +void init_basis_heading_and_non_basic_columns_vector(vector & basis, + vector & basis_heading, + vector & non_basic_columns) { + init_basic_part_of_basis_heading(basis, basis_heading); + init_non_basic_part_of_basis_heading(basis_heading, non_basic_columns); +} + +void change_basis(unsigned entering, unsigned leaving, vector& basis, vector& nbasis, vector & basis_heading) { + int place_in_basis = basis_heading[leaving]; + int place_in_non_basis = - basis_heading[entering] - 1; + basis_heading[entering] = place_in_basis; + basis_heading[leaving] = -place_in_non_basis - 1; + basis[place_in_basis] = entering; + nbasis[place_in_non_basis] = leaving; +} + + +#ifdef LEAN_DEBUG +void test_small_lu(lp_settings & settings) { + std::cout << " test_small_lu" << std::endl; + static_matrix m(3, 6); + vector basis(3); + basis[0] = 0; + basis[1] = 1; + basis[2] = 3; + + m(0, 0) = 1; m(0, 2)= 3.9; m(2, 3) = 11; m(0, 5) = -3; + m(1, 1) = 4; m(1, 4) = 7; + m(2, 0) = 1.8; m(2, 2) = 5; m(2, 4) = 2; m(2, 5) = 8; + +#ifdef LEAN_DEBUG + print_matrix(m, std::cout); +#endif + vector heading = allocate_basis_heading(m.column_count()); + vector non_basic_columns; + init_basis_heading_and_non_basic_columns_vector(basis, heading, non_basic_columns); + lu l(m, basis, settings); + lean_assert(l.is_correct(basis)); + indexed_vector w(m.row_count()); + std::cout << "entering 2, leaving 0" << std::endl; + l.prepare_entering(2, w); // to init vector w + l.replace_column(0, w, heading[0]); + change_basis(2, 0, basis, non_basic_columns, heading); + // #ifdef LEAN_DEBUG + // std::cout << "we were factoring " << std::endl; + // print_matrix(get_B(l)); + // #endif + lean_assert(l.is_correct(basis)); + std::cout << "entering 4, leaving 3" << std::endl; + l.prepare_entering(4, w); // to init vector w + l.replace_column(0, w, heading[3]); + change_basis(4, 3, basis, non_basic_columns, heading); + std::cout << "we were factoring " << std::endl; +#ifdef LEAN_DEBUG + { + auto bl = get_B(l, basis); + print_matrix(&bl, std::cout); + } +#endif + lean_assert(l.is_correct(basis)); + + std::cout << "entering 5, leaving 1" << std::endl; + l.prepare_entering(5, w); // to init vector w + l.replace_column(0, w, heading[1]); + change_basis(5, 1, basis, non_basic_columns, heading); + std::cout << "we were factoring " << std::endl; +#ifdef LEAN_DEBUG + { + auto bl = get_B(l, basis); + print_matrix(&bl, std::cout); + } +#endif + lean_assert(l.is_correct(basis)); + std::cout << "entering 3, leaving 2" << std::endl; + l.prepare_entering(3, w); // to init vector w + l.replace_column(0, w, heading[2]); + change_basis(3, 2, basis, non_basic_columns, heading); + std::cout << "we were factoring " << std::endl; +#ifdef LEAN_DEBUG + { + auto bl = get_B(l, basis); + print_matrix(&bl, std::cout); + } +#endif + lean_assert(l.is_correct(basis)); + + m.add_row(); + m.add_column(); + m.add_row(); + m.add_column(); + for (unsigned i = 0; i < m.column_count(); i++) { + m(3, i) = i; + m(4, i) = i * i; // to make the rows linearly independent + } + unsigned j = m.column_count() ; + basis.push_back(j-2); + heading.push_back(basis.size() - 1); + basis.push_back(j-1); + heading.push_back(basis.size() - 1); + + auto columns_to_replace = l.get_set_of_columns_to_replace_for_add_last_rows(heading); + l.add_last_rows_to_B(heading, columns_to_replace); + std::cout << "here" << std::endl; + lean_assert(l.is_correct(basis)); +} + +#endif + +void fill_long_row(sparse_matrix &m, int i) { + int n = m.dimension(); + for (int j = 0; j < n; j ++) { + m (i, (j + i) % n) = j * j; + } +} + +void fill_long_row(static_matrix &m, int i) { + int n = m.column_count(); + for (int j = 0; j < n; j ++) { + m (i, (j + i) % n) = j * j; + } +} + + +void fill_long_row_exp(sparse_matrix &m, int i) { + int n = m.dimension(); + + for (int j = 0; j < n; j ++) { + m(i, j) = my_random() % 20; + } +} + +void fill_long_row_exp(static_matrix &m, int i) { + int n = m.column_count(); + + for (int j = 0; j < n; j ++) { + m(i, j) = my_random() % 20; + } +} + +void fill_larger_sparse_matrix_exp(sparse_matrix & m){ + for ( unsigned i = 0; i < m.dimension(); i++ ) + fill_long_row_exp(m, i); +} + +void fill_larger_sparse_matrix_exp(static_matrix & m){ + for ( unsigned i = 0; i < m.row_count(); i++ ) + fill_long_row_exp(m, i); +} + + +void fill_larger_sparse_matrix(sparse_matrix & m){ + for ( unsigned i = 0; i < m.dimension(); i++ ) + fill_long_row(m, i); +} + +void fill_larger_sparse_matrix(static_matrix & m){ + for ( unsigned i = 0; i < m.row_count(); i++ ) + fill_long_row(m, i); +} + + +int perm_id = 0; + +#ifdef LEAN_DEBUG +void test_larger_lu_exp(lp_settings & settings) { + std::cout << " test_larger_lu_exp" << std::endl; + static_matrix m(6, 12); + vector basis(6); + basis[0] = 1; + basis[1] = 3; + basis[2] = 0; + basis[3] = 4; + basis[4] = 5; + basis[5] = 6; + + + fill_larger_sparse_matrix_exp(m); + // print_matrix(m); + vector heading = allocate_basis_heading(m.column_count()); + vector non_basic_columns; + init_basis_heading_and_non_basic_columns_vector(basis, heading, non_basic_columns); + lu l(m, basis, settings); + + dense_matrix left_side = l.get_left_side(basis); + dense_matrix right_side = l.get_right_side(); + lean_assert(left_side == right_side); + int leaving = 3; + int entering = 8; + for (unsigned i = 0; i < m.row_count(); i++) { + std::cout << static_cast(m(i, entering)) << std::endl; + } + + indexed_vector w(m.row_count()); + + l.prepare_entering(entering, w); + l.replace_column(0, w, heading[leaving]); + change_basis(entering, leaving, basis, non_basic_columns, heading); + lean_assert(l.is_correct(basis)); + + l.prepare_entering(11, w); // to init vector w + l.replace_column(0, w, heading[0]); + change_basis(11, 0, basis, non_basic_columns, heading); + lean_assert(l.is_correct(basis)); +} + +void test_larger_lu_with_holes(lp_settings & settings) { + std::cout << " test_larger_lu_with_holes" << std::endl; + static_matrix m(8, 9); + vector basis(8); + for (unsigned i = 0; i < m.row_count(); i++) { + basis[i] = i; + } + m(0, 0) = 1; m(0, 1) = 2; m(0, 2) = 3; m(0, 3) = 4; m(0, 4) = 5; m(0, 8) = 99; + /* */ m(1, 1) =- 6; m(1, 2) = 7; m(1, 3) = 8; m(1, 4) = 9; + /* */ m(2, 2) = 10; + /* */ m(3, 2) = 11; m(3, 3) = -12; + /* */ m(4, 2) = 13; m(4, 3) = 14; m(4, 4) = 15; + // the rest of the matrix is denser + m(5, 4) = 28; m(5, 5) = -18; m(5, 6) = 19; m(5, 7) = 25; + /* */ m(6, 5) = 20; m(6, 6) = -21; + /* */ m(7, 5) = 22; m(7, 6) = 23; m(7, 7) = 24; m(7, 8) = 88; + print_matrix(m, std::cout); + vector heading = allocate_basis_heading(m.column_count()); + vector non_basic_columns; + init_basis_heading_and_non_basic_columns_vector(basis, heading, non_basic_columns); + lu l(m, basis, settings); + std::cout << "printing factorization" << std::endl; + for (int i = l.tail_size() - 1; i >=0; i--) { + auto lp = l.get_lp_matrix(i); + lp->set_number_of_columns(m.row_count()); + lp->set_number_of_rows(m.row_count()); + print_matrix( lp, std::cout); + } + + dense_matrix left_side = l.get_left_side(basis); + dense_matrix right_side = l.get_right_side(); + if (!(left_side == right_side)) { + std::cout << "different sides" << std::endl; + } + + indexed_vector w(m.row_count()); + l.prepare_entering(8, w); // to init vector w + l.replace_column(0, w, heading[0]); + change_basis(8, 0, basis, non_basic_columns, heading); + lean_assert(l.is_correct(basis)); +} + + +void test_larger_lu(lp_settings& settings) { + std::cout << " test_larger_lu" << std::endl; + static_matrix m(6, 12); + vector basis(6); + basis[0] = 1; + basis[1] = 3; + basis[2] = 0; + basis[3] = 4; + basis[4] = 5; + basis[5] = 6; + + + fill_larger_sparse_matrix(m); + print_matrix(m, std::cout); + + vector heading = allocate_basis_heading(m.column_count()); + vector non_basic_columns; + init_basis_heading_and_non_basic_columns_vector(basis, heading, non_basic_columns); + auto l = lu (m, basis, settings); + // std::cout << "printing factorization" << std::endl; + // for (int i = lu.tail_size() - 1; i >=0; i--) { + // auto lp = lu.get_lp_matrix(i); + // lp->set_number_of_columns(m.row_count()); + // lp->set_number_of_rows(m.row_count()); + // print_matrix(* lp); + // } + + dense_matrix left_side = l.get_left_side(basis); + dense_matrix right_side = l.get_right_side(); + if (!(left_side == right_side)) { + std::cout << "left side" << std::endl; + print_matrix(&left_side, std::cout); + std::cout << "right side" << std::endl; + print_matrix(&right_side, std::cout); + + std::cout << "different sides" << std::endl; + std::cout << "initial factorization is incorrect" << std::endl; + exit(1); + } + indexed_vector w(m.row_count()); + l.prepare_entering(9, w); // to init vector w + l.replace_column(0, w, heading[0]); + change_basis(9, 0, basis, non_basic_columns, heading); + lean_assert(l.is_correct(basis)); +} + + +void test_lu(lp_settings & settings) { + test_small_lu(settings); + test_larger_lu(settings); + test_larger_lu_with_holes(settings); + test_larger_lu_exp(settings); +} +#endif + + + + + + +void init_b(vector & b, sparse_matrix & m, vector& x) { + for (unsigned i = 0; i < m.dimension(); i++) { + b.push_back(m.dot_product_with_row(i, x)); + } +} + +void init_b(vector & b, static_matrix & m, vector & x) { + for (unsigned i = 0; i < m.row_count(); i++) { + b.push_back(m.dot_product_with_row(i, x)); + } +} + + +void test_lp_0() { + std::cout << " test_lp_0 " << std::endl; + static_matrix m_(3, 7); + m_(0, 0) = 3; m_(0, 1) = 2; m_(0, 2) = 1; m_(0, 3) = 2; m_(0, 4) = 1; + m_(1, 0) = 1; m_(1, 1) = 1; m_(1, 2) = 1; m_(1, 3) = 1; m_(1, 5) = 1; + m_(2, 0) = 4; m_(2, 1) = 3; m_(2, 2) = 3; m_(2, 3) = 4; m_(2, 6) = 1; + vector x_star(7); + x_star[0] = 225; x_star[1] = 117; x_star[2] = 420; + x_star[3] = x_star[4] = x_star[5] = x_star[6] = 0; + vector b; + init_b(b, m_, x_star); + vector basis(3); + basis[0] = 0; basis[1] = 1; basis[2] = 2; + vector costs(7); + costs[0] = 19; + costs[1] = 13; + costs[2] = 12; + costs[3] = 17; + costs[4] = 0; + costs[5] = 0; + costs[6] = 0; + + vector column_types(7, column_type::low_bound); + vector upper_bound_values; + lp_settings settings; + simple_column_namer cn; + vector nbasis; + vector heading; + + lp_primal_core_solver lpsolver(m_, b, x_star, basis, nbasis, heading, costs, column_types, upper_bound_values, settings, cn); + + lpsolver.solve(); +} + +void test_lp_1() { + std::cout << " test_lp_1 " << std::endl; + static_matrix m(4, 7); + m(0, 0) = 1; m(0, 1) = 3; m(0, 2) = 1; m(0, 3) = 1; + m(1, 0) = -1; m(1, 2) = 3; m(1, 4) = 1; + m(2, 0) = 2; m(2, 1) = -1; m(2, 2) = 2; m(2, 5) = 1; + m(3, 0) = 2; m(3, 1) = 3; m(3, 2) = -1; m(3, 6) = 1; +#ifdef LEAN_DEBUG + print_matrix(m, std::cout); +#endif + vector x_star(7); + x_star[0] = 0; x_star[1] = 0; x_star[2] = 0; + x_star[3] = 3; x_star[4] = 2; x_star[5] = 4; x_star[6] = 2; + + vector basis(4); + basis[0] = 3; basis[1] = 4; basis[2] = 5; basis[3] = 6; + + vector b; + b.push_back(3); + b.push_back(2); + b.push_back(4); + b.push_back(2); + + vector costs(7); + costs[0] = 5; + costs[1] = 5; + costs[2] = 3; + costs[3] = 0; + costs[4] = 0; + costs[5] = 0; + costs[6] = 0; + + + + vector column_types(7, column_type::low_bound); + vector upper_bound_values; + + std::cout << "calling lp\n"; + lp_settings settings; + simple_column_namer cn; + + vector nbasis; + vector heading; + + lp_primal_core_solver lpsolver(m, b, + x_star, + basis, + nbasis, heading, + costs, + column_types, upper_bound_values, settings, cn); + + lpsolver.solve(); +} + + +void test_lp_primal_core_solver() { + test_lp_0(); + test_lp_1(); +} + + +#ifdef LEAN_DEBUG +template +void test_swap_rows_with_permutation(sparse_matrix& m){ + std::cout << "testing swaps" << std::endl; + unsigned dim = m.row_count(); + dense_matrix original(&m); + permutation_matrix q(dim); + print_matrix(m, std::cout); + lean_assert(original == q * m); + for (int i = 0; i < 100; i++) { + unsigned row1 = my_random() % dim; + unsigned row2 = my_random() % dim; + if (row1 == row2) continue; + std::cout << "swap " << row1 << " " << row2 << std::endl; + m.swap_rows(row1, row2); + q.transpose_from_left(row1, row2); + lean_assert(original == q * m); + print_matrix(m, std::cout); + std::cout << std::endl; + } +} +#endif +template +void fill_matrix(sparse_matrix& m); // forward definition +#ifdef LEAN_DEBUG +template +void test_swap_cols_with_permutation(sparse_matrix& m){ + std::cout << "testing swaps" << std::endl; + unsigned dim = m.row_count(); + dense_matrix original(&m); + permutation_matrix q(dim); + print_matrix(m, std::cout); + lean_assert(original == q * m); + for (int i = 0; i < 100; i++) { + unsigned row1 = my_random() % dim; + unsigned row2 = my_random() % dim; + if (row1 == row2) continue; + std::cout << "swap " << row1 << " " << row2 << std::endl; + m.swap_rows(row1, row2); + q.transpose_from_right(row1, row2); + lean_assert(original == q * m); + print_matrix(m, std::cout); + std::cout << std::endl; + } +} + + +template +void test_swap_rows(sparse_matrix& m, unsigned i0, unsigned i1){ + std::cout << "test_swap_rows(" << i0 << "," << i1 << ")" << std::endl; + sparse_matrix mcopy(m.dimension()); + for (unsigned i = 0; i < m.dimension(); i++) + for (unsigned j = 0; j < m.dimension(); j++) { + mcopy(i, j)= m(i, j); + } + std::cout << "swapping rows "<< i0 << "," << i1 << std::endl; + m.swap_rows(i0, i1); + + for (unsigned j = 0; j < m.dimension(); j++) { + lean_assert(mcopy(i0, j) == m(i1, j)); + lean_assert(mcopy(i1, j) == m(i0, j)); + } +} +template +void test_swap_columns(sparse_matrix& m, unsigned i0, unsigned i1){ + std::cout << "test_swap_columns(" << i0 << "," << i1 << ")" << std::endl; + sparse_matrix mcopy(m.dimension()); + for (unsigned i = 0; i < m.dimension(); i++) + for (unsigned j = 0; j < m.dimension(); j++) { + mcopy(i, j)= m(i, j); + } + m.swap_columns(i0, i1); + + for (unsigned j = 0; j < m.dimension(); j++) { + lean_assert(mcopy(j, i0) == m(j, i1)); + lean_assert(mcopy(j, i1) == m(j, i0)); + } + + for (unsigned i = 0; i < m.dimension(); i++) { + if (i == i0 || i == i1) + continue; + for (unsigned j = 0; j < m.dimension(); j++) { + lean_assert(mcopy(j, i)== m(j, i)); + } + } +} +#endif + +template +void fill_matrix(sparse_matrix& m){ + int v = 0; + for (int i = m.dimension() - 1; i >= 0; i--) { + for (int j = m.dimension() - 1; j >=0; j--){ + m(i, j) = v++; + } + } +} + +void test_pivot_like_swaps_and_pivot(){ + sparse_matrix m(10); + fill_matrix(m); + // print_matrix(m); +// pivot at 2,7 + m.swap_columns(0, 7); + // print_matrix(m); + m.swap_rows(2, 0); + // print_matrix(m); + for (unsigned i = 1; i < m.dimension(); i++) { + m(i, 0) = 0; + } + // print_matrix(m); + +// say pivot at 3,4 + m.swap_columns(1, 4); + // print_matrix(m); + m.swap_rows(1, 3); + // print_matrix(m); + + vector row; + double alpha = 2.33; + unsigned pivot_row = 1; + unsigned target_row = 2; + unsigned pivot_row_0 = 3; + double beta = 3.1; + m(target_row, 3) = 0; + m(target_row, 5) = 0; + m(pivot_row, 6) = 0; +#ifdef LEAN_DEBUG + print_matrix(m, std::cout); +#endif + + for (unsigned j = 0; j < m.dimension(); j++) { + row.push_back(m(target_row, j) + alpha * m(pivot_row, j) + beta * m(pivot_row_0, j)); + } + + for (auto & t : row) { + std::cout << t << ","; + } + std::cout << std::endl; + lp_settings settings; + m.pivot_row_to_row(pivot_row, alpha, target_row, settings); + m.pivot_row_to_row(pivot_row_0, beta, target_row, settings); + // print_matrix(m); + for (unsigned j = 0; j < m.dimension(); j++) { + lean_assert(abs(row[j] - m(target_row, j)) < 0.00000001); + } +} + +#ifdef LEAN_DEBUG +void test_swap_rows() { + sparse_matrix m(10); + fill_matrix(m); + // print_matrix(m); + test_swap_rows(m, 3, 5); + + test_swap_rows(m, 1, 3); + + + test_swap_rows(m, 1, 3); + + test_swap_rows(m, 1, 7); + + test_swap_rows(m, 3, 7); + + test_swap_rows(m, 0, 7); + + m(0, 4) = 1; + // print_matrix(m); + test_swap_rows(m, 0, 7); + +// go over some corner cases + sparse_matrix m0(2); + test_swap_rows(m0, 0, 1); + m0(0, 0) = 3; + test_swap_rows(m0, 0, 1); + m0(1, 0) = 3; + test_swap_rows(m0, 0, 1); + + + sparse_matrix m1(10); + test_swap_rows(m1, 0, 1); + m1(0, 0) = 3; + test_swap_rows(m1, 0, 1); + m1(1, 0) = 3; + m1(0, 3) = 5; + m1(1, 3) = 4; + m1(1, 8) = 8; + m1(1, 9) = 8; + + test_swap_rows(m1, 0, 1); + + sparse_matrix m2(3); + test_swap_rows(m2, 0, 1); + m2(0, 0) = 3; + test_swap_rows(m2, 0, 1); + m2(2, 0) = 3; + test_swap_rows(m2, 0, 2); +} + +void fill_uniformly(sparse_matrix & m, unsigned dim) { + int v = 0; + for (unsigned i = 0; i < dim; i++) { + for (unsigned j = 0; j < dim; j++) { + m(i, j) = v++; + } + } +} + +void fill_uniformly(dense_matrix & m, unsigned dim) { + int v = 0; + for (unsigned i = 0; i < dim; i++) { + for (unsigned j = 0; j < dim; j++) { + m.set_elem(i, j, v++); + } + } +} + +void sparse_matrix_with_permutaions_test() { + unsigned dim = 4; + sparse_matrix m(dim); + fill_uniformly(m, dim); + dense_matrix dm(dim, dim); + fill_uniformly(dm, dim); + dense_matrix dm0(dim, dim); + fill_uniformly(dm0, dim); + permutation_matrix q0(dim); + q0[0] = 1; + q0[1] = 0; + q0[2] = 3; + q0[3] = 2; + permutation_matrix q1(dim); + q1[0] = 1; + q1[1] = 2; + q1[2] = 3; + q1[3] = 0; + permutation_matrix p0(dim); + p0[0] = 1; + p0[1] = 0; + p0[2] = 3; + p0[3] = 2; + permutation_matrix p1(dim); + p1[0] = 1; + p1[1] = 2; + p1[2] = 3; + p1[3] = 0; + + m.multiply_from_left(q0); + for (unsigned i = 0; i < dim; i++) { + for (unsigned j = 0; j < dim; j++) { + lean_assert(m(i, j) == dm0.get_elem(q0[i], j)); + } + } + + auto q0_dm = q0 * dm; + lean_assert(m == q0_dm); + + m.multiply_from_left(q1); + for (unsigned i = 0; i < dim; i++) { + for (unsigned j = 0; j < dim; j++) { + lean_assert(m(i, j) == dm0.get_elem(q0[q1[i]], j)); + } + } + + + auto q1_q0_dm = q1 * q0_dm; + + lean_assert(m == q1_q0_dm); + + m.multiply_from_right(p0); + + for (unsigned i = 0; i < dim; i++) { + for (unsigned j = 0; j < dim; j++) { + lean_assert(m(i, j) == dm0.get_elem(q0[q1[i]], p0[j])); + } + } + + auto q1_q0_dm_p0 = q1_q0_dm * p0; + + lean_assert(m == q1_q0_dm_p0); + + m.multiply_from_right(p1); + + for (unsigned i = 0; i < dim; i++) { + for (unsigned j = 0; j < dim; j++) { + lean_assert(m(i, j) == dm0.get_elem(q0[q1[i]], p1[p0[j]])); + } + } + + auto q1_q0_dm_p0_p1 = q1_q0_dm_p0 * p1; + lean_assert(m == q1_q0_dm_p0_p1); + + m.multiply_from_right(p1); + for (unsigned i = 0; i < dim; i++) { + for (unsigned j = 0; j < dim; j++) { + lean_assert(m(i, j) == dm0.get_elem(q0[q1[i]], p1[p1[p0[j]]])); + } + } + auto q1_q0_dm_p0_p1_p1 = q1_q0_dm_p0_p1 * p1; + + lean_assert(m == q1_q0_dm_p0_p1_p1); +} + +void test_swap_columns() { + sparse_matrix m(10); + fill_matrix(m); + // print_matrix(m); + + test_swap_columns(m, 3, 5); + + test_swap_columns(m, 1, 3); + + test_swap_columns(m, 1, 3); + + // print_matrix(m); + test_swap_columns(m, 1, 7); + + test_swap_columns(m, 3, 7); + + test_swap_columns(m, 0, 7); + + test_swap_columns(m, 0, 7); + +// go over some corner cases + sparse_matrix m0(2); + test_swap_columns(m0, 0, 1); + m0(0, 0) = 3; + test_swap_columns(m0, 0, 1); + m0(0, 1) = 3; + test_swap_columns(m0, 0, 1); + + + sparse_matrix m1(10); + test_swap_columns(m1, 0, 1); + m1(0, 0) = 3; + test_swap_columns(m1, 0, 1); + m1(0, 1) = 3; + m1(3, 0) = 5; + m1(3, 1) = 4; + m1(8, 1) = 8; + m1(9, 1) = 8; + + test_swap_columns(m1, 0, 1); + + sparse_matrix m2(3); + test_swap_columns(m2, 0, 1); + m2(0, 0) = 3; + test_swap_columns(m2, 0, 1); + m2(0, 2) = 3; + test_swap_columns(m2, 0, 2); +} + + + +void test_swap_operations() { + test_swap_rows(); + test_swap_columns(); +} + +void test_dense_matrix() { + dense_matrix d(3, 2); + d.set_elem(0, 0, 1); + d.set_elem(1, 1, 2); + d.set_elem(2, 0, 3); + // print_matrix(d); + + dense_matrix unit(2, 2); + d.set_elem(0, 0, 1); + d.set_elem(1, 1, 1); + + dense_matrix c = d * unit; + + // print_matrix(d); + + dense_matrix perm(3, 3); + perm.set_elem(0, 1, 1); + perm.set_elem(1, 0, 1); + perm.set_elem(2, 2, 1); + auto c1 = perm * d; + // print_matrix(c1); + + + dense_matrix p2(2, 2); + p2.set_elem(0, 1, 1); + p2.set_elem(1, 0, 1); + auto c2 = d * p2; +} +#endif + + + +vector> vector_of_permutaions() { + vector> ret; + { + permutation_matrix p0(5); + p0[0] = 1; p0[1] = 2; p0[2] = 3; p0[3] = 4; + p0[4] = 0; + ret.push_back(p0); + } + { + permutation_matrix p0(5); + p0[0] = 2; p0[1] = 0; p0[2] = 1; p0[3] = 4; + p0[4] = 3; + ret.push_back(p0); + } + return ret; +} + +void test_apply_reverse_from_right_to_perm(permutation_matrix & l) { + permutation_matrix p(5); + p[0] = 4; p[1] = 2; p[2] = 0; p[3] = 3; + p[4] = 1; + + permutation_matrix pclone(5); + pclone[0] = 4; pclone[1] = 2; pclone[2] = 0; pclone[3] = 3; + pclone[4] = 1; + + p.multiply_by_reverse_from_right(l); +#ifdef LEAN_DEBUG + auto rev = l.get_inverse(); + auto rs = pclone * rev; + lean_assert(p == rs) +#endif +} + +void test_apply_reverse_from_right() { + auto vec = vector_of_permutaions(); + for (unsigned i = 0; i < vec.size(); i++) { + test_apply_reverse_from_right_to_perm(vec[i]); + } +} + +void test_permutations() { + std::cout << "test permutations" << std::endl; + test_apply_reverse_from_right(); + vector v; v.resize(5, 0); + v[1] = 1; + v[3] = 3; + permutation_matrix p(5); + p[0] = 4; p[1] = 2; p[2] = 0; p[3] = 3; + p[4] = 1; + + indexed_vector vi(5); + vi.set_value(1, 1); + vi.set_value(3, 3); + + p.apply_reverse_from_right_to_T(v); + p.apply_reverse_from_right_to_T(vi); + lean_assert(vectors_are_equal(v, vi.m_data)); + lean_assert(vi.is_OK()); +} + +void lp_solver_test() { + // lp_revised_solver lp_revised; + // lp_revised.get_minimal_solution(); +} + +bool get_int_from_args_parser(const char * option, argument_parser & args_parser, unsigned & n) { + std::string s = args_parser.get_option_value(option); + if (s.size() > 0) { + n = atoi(s.c_str()); + return true; + } + return false; +} + +bool get_double_from_args_parser(const char * option, argument_parser & args_parser, double & n) { + std::string s = args_parser.get_option_value(option); + if (s.size() > 0) { + n = atof(s.c_str()); + return true; + } + return false; +} + + +void update_settings(argument_parser & args_parser, lp_settings& settings) { + unsigned n; + settings.m_simplex_strategy = simplex_strategy_enum::no_tableau; + if (get_int_from_args_parser("--rep_frq", args_parser, n)) + settings.report_frequency = n; + else + settings.report_frequency = args_parser.option_is_used("--mpq")? 80: 1000; + + settings.print_statistics = true; + + if (get_int_from_args_parser("--percent_for_enter", args_parser, n)) + settings.percent_of_entering_to_check = n; + if (get_int_from_args_parser("--partial_pivot", args_parser, n)) { + std::cout << "setting partial pivot constant to " << n << std::endl; + settings.c_partial_pivoting = n; + } + if (get_int_from_args_parser("--density", args_parser, n)) { + double density = static_cast(n) / 100.0; + std::cout << "setting density to " << density << std::endl; + settings.density_threshold = density; + } + if (get_int_from_args_parser("--maxng", args_parser, n)) + settings.max_number_of_iterations_with_no_improvements = n; + double d; + if (get_double_from_args_parser("--harris_toler", args_parser, d)) { + std::cout << "setting harris_feasibility_tolerance to " << d << std::endl; + settings.harris_feasibility_tolerance = d; + } + if (get_int_from_args_parser("--random_seed", args_parser, n)) { + settings.random_seed = n; + } + if (get_int_from_args_parser("--simplex_strategy", args_parser, n)) { + settings.simplex_strategy() = static_cast(n); + } +} + +template +void setup_solver(unsigned max_iterations, unsigned time_limit, bool look_for_min, argument_parser & args_parser, lp_solver * solver) { + if (max_iterations > 0) + solver->set_max_iterations_per_stage(max_iterations); + + if (time_limit > 0) + solver->set_time_limit(time_limit); + + if (look_for_min) + solver->flip_costs(); + + update_settings(args_parser, solver->settings()); +} + +bool values_are_one_percent_close(double a, double b); + +void print_x(mps_reader & reader, lp_solver * solver) { + for (auto name : reader.column_names()) { + std::cout << name << "=" << solver->get_column_value_by_name(name) << ' '; + } + std::cout << std::endl; +} + +void compare_solutions(mps_reader & reader, lp_solver * solver, lp_solver * solver0) { + for (auto name : reader.column_names()) { + double a = solver->get_column_value_by_name(name); + double b = solver0->get_column_value_by_name(name); + if (!values_are_one_percent_close(a, b)) { + std::cout << "different values for " << name << ":" << a << " and " << b << std::endl; + } + } +} + + +void solve_mps_double(std::string file_name, bool look_for_min, unsigned max_iterations, unsigned time_limit, bool dual, bool compare_with_primal, argument_parser & args_parser) { + mps_reader reader(file_name); + reader.read(); + if (!reader.is_ok()) { + std::cout << "cannot process " << file_name << std::endl; + return; + } + + lp_solver * solver = reader.create_solver(dual); + setup_solver(max_iterations, time_limit, look_for_min, args_parser, solver); + int begin = get_millisecond_count(); + if (dual) { + std::cout << "solving for dual" << std::endl; + } + solver->find_maximal_solution(); + int span = get_millisecond_span(begin); + std::cout << "Status: " << lp_status_to_string(solver->get_status()) << std::endl; + if (solver->get_status() == lp_status::OPTIMAL) { + if (reader.column_names().size() < 20) { + print_x(reader, solver); + } + double cost = solver->get_current_cost(); + if (look_for_min) { + cost = -cost; + } + std::cout << "cost = " << cost << std::endl; + } + std::cout << "processed in " << span / 1000.0 << " seconds, running for " << solver->m_total_iterations << " iterations" << " one iter for " << (double)span/solver->m_total_iterations << " ms" << std::endl; + if (compare_with_primal) { + auto * primal_solver = reader.create_solver(false); + setup_solver(max_iterations, time_limit, look_for_min, args_parser, primal_solver); + primal_solver->find_maximal_solution(); + if (solver->get_status() != primal_solver->get_status()) { + std::cout << "statuses are different: dual " << lp_status_to_string(solver->get_status()) << " primal = " << lp_status_to_string(primal_solver->get_status()) << std::endl; + } else { + if (solver->get_status() == lp_status::OPTIMAL) { + double cost = solver->get_current_cost(); + if (look_for_min) { + cost = -cost; + } + double primal_cost = primal_solver->get_current_cost(); + if (look_for_min) { + primal_cost = -primal_cost; + } + std::cout << "primal cost = " << primal_cost << std::endl; + if (!values_are_one_percent_close(cost, primal_cost)) { + compare_solutions(reader, primal_solver, solver); + print_x(reader, primal_solver); + std::cout << "dual cost is " << cost << ", but primal cost is " << primal_cost << std::endl; + lean_assert(false); + } + } + } + delete primal_solver; + } + delete solver; +} + +void solve_mps_rational(std::string file_name, bool look_for_min, unsigned max_iterations, unsigned time_limit, bool dual, argument_parser & args_parser) { + mps_reader reader(file_name); + reader.read(); + if (reader.is_ok()) { + auto * solver = reader.create_solver(dual); + setup_solver(max_iterations, time_limit, look_for_min, args_parser, solver); + int begin = get_millisecond_count(); + solver->find_maximal_solution(); + std::cout << "Status: " << lp_status_to_string(solver->get_status()) << std::endl; + + if (solver->get_status() == lp_status::OPTIMAL) { + // for (auto name: reader.column_names()) { + // std::cout << name << "=" << solver->get_column_value_by_name(name) << ' '; + // } + lean::mpq cost = solver->get_current_cost(); + if (look_for_min) { + cost = -cost; + } + std::cout << "cost = " << cost.get_double() << std::endl; + } + std::cout << "processed in " << get_millisecond_span(begin) / 1000.0 << " seconds, running for " << solver->m_total_iterations << " iterations" << std::endl; + delete solver; + } else { + std::cout << "cannot process " << file_name << std::endl; + } +} +void get_time_limit_and_max_iters_from_parser(argument_parser & args_parser, unsigned & time_limit, unsigned & max_iters); // forward definition + +void solve_mps(std::string file_name, bool look_for_min, unsigned max_iterations, unsigned time_limit, bool solve_for_rational, bool dual, bool compare_with_primal, argument_parser & args_parser) { + if (!solve_for_rational) { + std::cout << "solving " << file_name << std::endl; + solve_mps_double(file_name, look_for_min, max_iterations, time_limit, dual, compare_with_primal, args_parser); + } + else { + std::cout << "solving " << file_name << " in rationals " << std::endl; + solve_mps_rational(file_name, look_for_min, max_iterations, time_limit, dual, args_parser); + } +} + +void solve_mps(std::string file_name, argument_parser & args_parser) { + bool look_for_min = args_parser.option_is_used("--min"); + unsigned max_iterations, time_limit; + bool solve_for_rational = args_parser.option_is_used("--mpq"); + bool dual = args_parser.option_is_used("--dual"); + bool compare_with_primal = args_parser.option_is_used("--compare_with_primal"); + get_time_limit_and_max_iters_from_parser(args_parser, time_limit, max_iterations); + solve_mps(file_name, look_for_min, max_iterations, time_limit, solve_for_rational, dual, compare_with_primal, args_parser); +} + +void solve_mps_in_rational(std::string file_name, bool dual, argument_parser & /*args_parser*/) { + std::cout << "solving " << file_name << std::endl; + + mps_reader reader(file_name); + reader.read(); + if (reader.is_ok()) { + auto * solver = reader.create_solver(dual); + solver->find_maximal_solution(); + std::cout << "status is " << lp_status_to_string(solver->get_status()) << std::endl; + if (solver->get_status() == lp_status::OPTIMAL) { + if (reader.column_names().size() < 20) { + for (auto name : reader.column_names()) { + std::cout << name << "=" << solver->get_column_value_by_name(name).get_double() << ' '; + } + } + std::cout << std::endl << "cost = " << numeric_traits::get_double(solver->get_current_cost()) << std::endl; + } + delete solver; + } else { + std::cout << "cannot process " << file_name << std::endl; + } +} + +void test_upair_queue() { + int n = 10; + binary_heap_upair_queue q(2); + std::unordered_map m; + for (int k = 0; k < 100; k++) { + int i = my_random()%n; + int j = my_random()%n; + q.enqueue(i, j, my_random()%n); + } + + q.remove(5, 5); + + while (!q.is_empty()) { + unsigned i, j; + q.dequeue(i, j); + } +} + +void test_binary_priority_queue() { + std::cout << "testing binary_heap_priority_queue..."; + auto q = binary_heap_priority_queue(10); + q.enqueue(2, 2); + q.enqueue(1, 1); + q.enqueue(9, 9); + q.enqueue(8, 8); + q.enqueue(5, 25); + q.enqueue(3, 3); + q.enqueue(4, 4); + q.enqueue(7, 30); + q.enqueue(6, 6); + q.enqueue(0, 0); + q.enqueue(5, 5); + q.enqueue(7, 7); + + for (unsigned i = 0; i < 10; i++) { + unsigned de = q.dequeue(); + lean_assert(i == de); + std::cout << de << std::endl; + } + q.enqueue(2, 2); + q.enqueue(1, 1); + q.enqueue(9, 9); + q.enqueue(8, 8); + q.enqueue(5, 5); + q.enqueue(3, 3); + q.enqueue(4, 4); + q.enqueue(7, 2); + q.enqueue(0, 1); + q.enqueue(6, 6); + q.enqueue(7, 7); + q.enqueue(33, 1000); + q.enqueue(20, 0); + q.dequeue(); + q.remove(33); + q.enqueue(0, 0); +#ifdef LEAN_DEBUG + unsigned t = 0; + while (q.size() > 0) { + unsigned d =q.dequeue(); + lean_assert(t++ == d); + std::cout << d << std::endl; + } +#endif + test_upair_queue(); + std::cout << " done" << std::endl; +} + +bool solution_is_feasible(std::string file_name, const std::unordered_map & solution) { + mps_reader reader(file_name); + reader.read(); + if (reader.is_ok()) { + lp_primal_simplex * solver = static_cast *>(reader.create_solver(false)); + return solver->solution_is_feasible(solution); + } + return false; +} + + +void solve_mps_with_known_solution(std::string file_name, std::unordered_map * solution, lp_status status, bool dual) { + std::cout << "solving " << file_name << std::endl; + mps_reader reader(file_name); + reader.read(); + if (reader.is_ok()) { + auto * solver = reader.create_solver(dual); + solver->find_maximal_solution(); + std::cout << "status is " << lp_status_to_string(solver->get_status()) << std::endl; + if (status != solver->get_status()){ + std::cout << "status should be " << lp_status_to_string(status) << std::endl; + lean_assert(status == solver->get_status()); + throw "status is wrong"; + } + if (solver->get_status() == lp_status::OPTIMAL) { + std::cout << "cost = " << solver->get_current_cost() << std::endl; + if (solution != nullptr) { + for (auto it : *solution) { + if (fabs(it.second - solver->get_column_value_by_name(it.first)) >= 0.000001) { + std::cout << "expected:" << it.first << "=" << + it.second <<", got " << solver->get_column_value_by_name(it.first) << std::endl; + } + lean_assert(fabs(it.second - solver->get_column_value_by_name(it.first)) < 0.000001); + } + } + if (reader.column_names().size() < 20) { + for (auto name : reader.column_names()) { + std::cout << name << "=" << solver->get_column_value_by_name(name) << ' '; + } + std::cout << std::endl; + } + } + delete solver; + } else { + std::cout << "cannot process " << file_name << std::endl; + } +} + +int get_random_rows() { + return 5 + my_random() % 2; +} + +int get_random_columns() { + return 5 + my_random() % 3; +} + +int get_random_int() { + return -1 + my_random() % 2; // (1.0 + RAND_MAX); +} + +void add_random_row(lp_primal_simplex * solver, int cols, int row) { + solver->add_constraint(lp_relation::Greater_or_equal, 1, row); + for (int i = 0; i < cols; i++) { + solver->set_row_column_coefficient(row, i, get_random_int()); + } +} + +void add_random_cost(lp_primal_simplex * solver, int cols) { + for (int i = 0; i < cols; i++) { + solver->set_cost_for_column(i, get_random_int()); + } +} + +lp_primal_simplex * generate_random_solver() { + int rows = get_random_rows(); + int cols = get_random_columns(); + auto * solver = new lp_primal_simplex(); + for (int i = 0; i < rows; i++) { + add_random_row(solver, cols, i); + } + add_random_cost(solver, cols); + return solver; +} + + + +void random_test_on_i(unsigned i) { + if (i % 1000 == 0) { + std::cout << "."; + } + srand(i); + auto *solver = generate_random_solver(); + solver->find_maximal_solution(); + // std::cout << lp_status_to_string(solver->get_status()) << std::endl; + delete solver; +} + +void random_test() { + for (unsigned i = 0; i < std::numeric_limits::max(); i++) { + try { + random_test_on_i(i); + } + catch (const char * error) { + std::cout << "i = " << i << ", throwing at ' " << error << "'" << std::endl; + break; + } + } +} + +#if _LINUX_ +void fill_file_names(vector &file_names, std::set & minimums) { + char *home_dir = getenv("HOME"); + if (home_dir == nullptr) { + std::cout << "cannot find home directory, don't know how to find the files"; + return; + } + std::string home_dir_str(home_dir); + file_names.push_back(home_dir_str + "/projects/lean/src/tests/util/lp/l0redund.mps"); + file_names.push_back(home_dir_str + "/projects/lean/src/tests/util/lp/l1.mps"); + file_names.push_back(home_dir_str + "/projects/lean/src/tests/util/lp/l2.mps"); + file_names.push_back(home_dir_str + "/projects/lean/src/tests/util/lp/l3.mps"); + file_names.push_back(home_dir_str + "/projects/lean/src/tests/util/lp/l4.mps"); + file_names.push_back(home_dir_str + "/projects/lean/src/tests/util/lp/l4fix.mps"); + file_names.push_back(home_dir_str + "/projects/lean/src/tests/util/lp/plan.mps"); + file_names.push_back(home_dir_str + "/projects/lean/src/tests/util/lp/samp2.mps"); + file_names.push_back(home_dir_str + "/projects/lean/src/tests/util/lp/murtagh.mps"); + file_names.push_back(home_dir_str + "/projects/lean/src/tests/util/lp/l0.mps"); + file_names.push_back(home_dir_str + "/projects/lean/src/tests/util/lp/test_files/netlib/AFIRO.SIF"); + file_names.push_back(home_dir_str + "/projects/lean/src/tests/util/lp/test_files/netlib/SC50B.SIF"); + file_names.push_back(home_dir_str + "/projects/lean/src/tests/util/lp/test_files/netlib/SC50A.SIF"); + file_names.push_back(home_dir_str + "/projects/lean/src/tests/util/lp/test_files/netlib/KB2.SIF"); + file_names.push_back(home_dir_str + "/projects/lean/src/tests/util/lp/test_files/netlib/SC105.SIF"); + file_names.push_back(home_dir_str + "/projects/lean/src/tests/util/lp/test_files/netlib/STOCFOR1.SIF"); + file_names.push_back(home_dir_str + "/projects/lean/src/tests/util/lp/test_files/netlib/ADLITTLE.SIF"); + file_names.push_back(home_dir_str + "/projects/lean/src/tests/util/lp/test_files/netlib/BLEND.SIF"); + file_names.push_back(home_dir_str + "/projects/lean/src/tests/util/lp/test_files/netlib/SCAGR7.SIF"); + file_names.push_back(home_dir_str + "/projects/lean/src/tests/util/lp/test_files/netlib/SC205.SIF"); + file_names.push_back(home_dir_str + "/projects/lean/src/tests/util/lp/test_files/netlib/SHARE2B.SIF"); + file_names.push_back(home_dir_str + "/projects/lean/src/tests/util/lp/test_files/netlib/RECIPELP.SIF"); + file_names.push_back(home_dir_str + "/projects/lean/src/tests/util/lp/test_files/netlib/LOTFI.SIF"); + file_names.push_back(home_dir_str + "/projects/lean/src/tests/util/lp/test_files/netlib/VTP-BASE.SIF"); + file_names.push_back(home_dir_str + "/projects/lean/src/tests/util/lp/test_files/netlib/SHARE1B.SIF"); + file_names.push_back(home_dir_str + "/projects/lean/src/tests/util/lp/test_files/netlib/BOEING2.SIF"); + file_names.push_back(home_dir_str + "/projects/lean/src/tests/util/lp/test_files/netlib/BORE3D.SIF"); + file_names.push_back(home_dir_str + "/projects/lean/src/tests/util/lp/test_files/netlib/SCORPION.SIF"); + file_names.push_back(home_dir_str + "/projects/lean/src/tests/util/lp/test_files/netlib/CAPRI.SIF"); + file_names.push_back(home_dir_str + "/projects/lean/src/tests/util/lp/test_files/netlib/BRANDY.SIF"); + file_names.push_back(home_dir_str + "/projects/lean/src/tests/util/lp/test_files/netlib/SCAGR25.SIF"); + file_names.push_back(home_dir_str + "/projects/lean/src/tests/util/lp/test_files/netlib/SCTAP1.SIF"); + file_names.push_back(home_dir_str + "/projects/lean/src/tests/util/lp/test_files/netlib/ISRAEL.SIF"); + file_names.push_back(home_dir_str + "/projects/lean/src/tests/util/lp/test_files/netlib/SCFXM1.SIF"); + file_names.push_back(home_dir_str + "/projects/lean/src/tests/util/lp/test_files/netlib/BANDM.SIF"); + file_names.push_back(home_dir_str + "/projects/lean/src/tests/util/lp/test_files/netlib/E226.SIF"); + file_names.push_back(home_dir_str + "/projects/lean/src/tests/util/lp/test_files/netlib/AGG.SIF"); + file_names.push_back(home_dir_str + "/projects/lean/src/tests/util/lp/test_files/netlib/GROW7.SIF"); + file_names.push_back(home_dir_str + "/projects/lean/src/tests/util/lp/test_files/netlib/ETAMACRO.SIF"); + file_names.push_back(home_dir_str + "/projects/lean/src/tests/util/lp/test_files/netlib/FINNIS.SIF"); + file_names.push_back(home_dir_str + "/projects/lean/src/tests/util/lp/test_files/netlib/SCSD1.SIF"); + file_names.push_back(home_dir_str + "/projects/lean/src/tests/util/lp/test_files/netlib/STANDATA.SIF"); + file_names.push_back(home_dir_str + "/projects/lean/src/tests/util/lp/test_files/netlib/STANDGUB.SIF"); + file_names.push_back(home_dir_str + "/projects/lean/src/tests/util/lp/test_files/netlib/BEACONFD.SIF"); + file_names.push_back(home_dir_str + "/projects/lean/src/tests/util/lp/test_files/netlib/STAIR.SIF"); + file_names.push_back(home_dir_str + "/projects/lean/src/tests/util/lp/test_files/netlib/STANDMPS.SIF"); + file_names.push_back(home_dir_str + "/projects/lean/src/tests/util/lp/test_files/netlib/GFRD-PNC.SIF"); + file_names.push_back(home_dir_str + "/projects/lean/src/tests/util/lp/test_files/netlib/SCRS8.SIF"); + file_names.push_back(home_dir_str + "/projects/lean/src/tests/util/lp/test_files/netlib/BOEING1.SIF"); + file_names.push_back(home_dir_str + "/projects/lean/src/tests/util/lp/test_files/netlib/MODSZK1.SIF"); + file_names.push_back(home_dir_str + "/projects/lean/src/tests/util/lp/test_files/netlib/DEGEN2.SIF"); + file_names.push_back(home_dir_str + "/projects/lean/src/tests/util/lp/test_files/netlib/FORPLAN.SIF"); + file_names.push_back(home_dir_str + "/projects/lean/src/tests/util/lp/test_files/netlib/AGG2.SIF"); + file_names.push_back(home_dir_str + "/projects/lean/src/tests/util/lp/test_files/netlib/AGG3.SIF"); + file_names.push_back(home_dir_str + "/projects/lean/src/tests/util/lp/test_files/netlib/SCFXM2.SIF"); + file_names.push_back(home_dir_str + "/projects/lean/src/tests/util/lp/test_files/netlib/SHELL.SIF"); + file_names.push_back(home_dir_str + "/projects/lean/src/tests/util/lp/test_files/netlib/PILOT4.SIF"); + file_names.push_back(home_dir_str + "/projects/lean/src/tests/util/lp/test_files/netlib/SCSD6.SIF"); + file_names.push_back(home_dir_str + "/projects/lean/src/tests/util/lp/test_files/netlib/SHIP04S.SIF"); + file_names.push_back(home_dir_str + "/projects/lean/src/tests/util/lp/test_files/netlib/SEBA.SIF"); + file_names.push_back(home_dir_str + "/projects/lean/src/tests/util/lp/test_files/netlib/GROW15.SIF"); + file_names.push_back(home_dir_str + "/projects/lean/src/tests/util/lp/test_files/netlib/FFFFF800.SIF"); + file_names.push_back(home_dir_str + "/projects/lean/src/tests/util/lp/test_files/netlib/BNL1.SIF"); + file_names.push_back(home_dir_str + "/projects/lean/src/tests/util/lp/test_files/netlib/PEROLD.SIF"); + file_names.push_back(home_dir_str + "/projects/lean/src/tests/util/lp/test_files/netlib/QAP8.SIF"); + file_names.push_back(home_dir_str + "/projects/lean/src/tests/util/lp/test_files/netlib/SCFXM3.SIF"); + file_names.push_back(home_dir_str + "/projects/lean/src/tests/util/lp/test_files/netlib/SHIP04L.SIF"); + file_names.push_back(home_dir_str + "/projects/lean/src/tests/util/lp/test_files/netlib/GANGES.SIF"); + file_names.push_back(home_dir_str + "/projects/lean/src/tests/util/lp/test_files/netlib/SCTAP2.SIF"); + file_names.push_back(home_dir_str + "/projects/lean/src/tests/util/lp/test_files/netlib/GROW22.SIF"); + file_names.push_back(home_dir_str + "/projects/lean/src/tests/util/lp/test_files/netlib/SHIP08S.SIF"); + file_names.push_back(home_dir_str + "/projects/lean/src/tests/util/lp/test_files/netlib/PILOT-WE.SIF"); + file_names.push_back(home_dir_str + "/projects/lean/src/tests/util/lp/test_files/netlib/MAROS.SIF"); + file_names.push_back(home_dir_str + "/projects/lean/src/tests/util/lp/test_files/netlib/STOCFOR2.SIF"); + file_names.push_back(home_dir_str + "/projects/lean/src/tests/util/lp/test_files/netlib/25FV47.SIF"); + file_names.push_back(home_dir_str + "/projects/lean/src/tests/util/lp/test_files/netlib/SHIP12S.SIF"); + file_names.push_back(home_dir_str + "/projects/lean/src/tests/util/lp/test_files/netlib/SCSD8.SIF"); + file_names.push_back(home_dir_str + "/projects/lean/src/tests/util/lp/test_files/netlib/FIT1P.SIF"); + file_names.push_back(home_dir_str + "/projects/lean/src/tests/util/lp/test_files/netlib/SCTAP3.SIF"); + file_names.push_back(home_dir_str + "/projects/lean/src/tests/util/lp/test_files/netlib/SIERRA.SIF"); + file_names.push_back(home_dir_str + "/projects/lean/src/tests/util/lp/test_files/netlib/PILOTNOV.SIF"); + file_names.push_back(home_dir_str + "/projects/lean/src/tests/util/lp/test_files/netlib/CZPROB.SIF"); + file_names.push_back(home_dir_str + "/projects/lean/src/tests/util/lp/test_files/netlib/FIT1D.SIF"); + file_names.push_back(home_dir_str + "/projects/lean/src/tests/util/lp/test_files/netlib/PILOT-JA.SIF"); + file_names.push_back(home_dir_str + "/projects/lean/src/tests/util/lp/test_files/netlib/SHIP08L.SIF"); + file_names.push_back(home_dir_str + "/projects/lean/src/tests/util/lp/test_files/netlib/BNL2.SIF"); + file_names.push_back(home_dir_str + "/projects/lean/src/tests/util/lp/test_files/netlib/NESM.SIF"); + file_names.push_back(home_dir_str + "/projects/lean/src/tests/util/lp/test_files/netlib/CYCLE.SIF"); + file_names.push_back(home_dir_str + "/projects/lean/src/tests/util/lp/test_files/acc-tight5.mps"); + file_names.push_back(home_dir_str + "/projects/lean/src/tests/util/lp/test_files/netlib/SHIP12L.SIF"); + file_names.push_back(home_dir_str + "/projects/lean/src/tests/util/lp/test_files/netlib/DEGEN3.SIF"); + file_names.push_back(home_dir_str + "/projects/lean/src/tests/util/lp/test_files/netlib/GREENBEA.SIF"); + file_names.push_back(home_dir_str + "/projects/lean/src/tests/util/lp/test_files/netlib/GREENBEB.SIF"); + file_names.push_back(home_dir_str + "/projects/lean/src/tests/util/lp/test_files/netlib/80BAU3B.SIF"); + file_names.push_back(home_dir_str + "/projects/lean/src/tests/util/lp/test_files/netlib/TRUSS.SIF"); + file_names.push_back(home_dir_str + "/projects/lean/src/tests/util/lp/test_files/netlib/D2Q06C.SIF"); + file_names.push_back(home_dir_str + "/projects/lean/src/tests/util/lp/test_files/netlib/WOODW.SIF"); + file_names.push_back(home_dir_str + "/projects/lean/src/tests/util/lp/test_files/netlib/QAP12.SIF"); + file_names.push_back(home_dir_str + "/projects/lean/src/tests/util/lp/test_files/netlib/D6CUBE.SIF"); + file_names.push_back(home_dir_str + "/projects/lean/src/tests/util/lp/test_files/netlib/PILOT.SIF"); + file_names.push_back(home_dir_str + "/projects/lean/src/tests/util/lp/test_files/netlib/DFL001.SIF"); + file_names.push_back(home_dir_str + "/projects/lean/src/tests/util/lp/test_files/netlib/WOOD1P.SIF"); + file_names.push_back(home_dir_str + "/projects/lean/src/tests/util/lp/test_files/netlib/FIT2P.SIF"); + file_names.push_back(home_dir_str + "/projects/lean/src/tests/util/lp/test_files/netlib/PILOT87.SIF"); + file_names.push_back(home_dir_str + "/projects/lean/src/tests/util/lp/test_files/netlib/STOCFOR3.SIF"); + file_names.push_back(home_dir_str + "/projects/lean/src/tests/util/lp/test_files/netlib/QAP15.SIF"); + file_names.push_back(home_dir_str + "/projects/lean/src/tests/util/lp/test_files/netlib/FIT2D.SIF"); + file_names.push_back(home_dir_str + "/projects/lean/src/tests/util/lp/test_files/netlib/MAROS-R7.SIF"); + minimums.insert("/projects/lean/src/tests/util/lp/test_files/netlib/FIT2P.SIF"); + minimums.insert("/projects/lean/src/tests/util/lp/test_files/netlib/DFL001.SIF"); + minimums.insert("/projects/lean/src/tests/util/lp/test_files/netlib/D2Q06C.SIF"); + minimums.insert("/projects/lean/src/tests/util/lp/test_files/netlib/80BAU3B.SIF"); + minimums.insert("/projects/lean/src/tests/util/lp/test_files/netlib/GREENBEB.SIF"); + minimums.insert("/projects/lean/src/tests/util/lp/test_files/netlib/GREENBEA.SIF"); + minimums.insert("/projects/lean/src/tests/util/lp/test_files/netlib/BNL2.SIF"); + minimums.insert("/projects/lean/src/tests/util/lp/test_files/netlib/SHIP08L.SIF"); + minimums.insert("/projects/lean/src/tests/util/lp/test_files/netlib/FIT1D.SIF"); + minimums.insert("/projects/lean/src/tests/util/lp/test_files/netlib/SCTAP3.SIF"); + minimums.insert("/projects/lean/src/tests/util/lp/test_files/netlib/SCSD8.SIF"); + minimums.insert("/projects/lean/src/tests/util/lp/test_files/netlib/SCSD6.SIF"); + minimums.insert("/projects/lean/src/tests/util/lp/test_files/netlib/MAROS-R7.SIF"); +} + +void test_out_dir(std::string out_dir) { + auto *out_dir_p = opendir(out_dir.c_str()); + if (out_dir_p == nullptr) { + std::cout << "creating directory " << out_dir << std::endl; +#ifdef LEAN_WINDOWS + int res = mkdir(out_dir.c_str()); +#else + int res = mkdir(out_dir.c_str(), S_IRWXU | S_IRWXG | S_IROTH | S_IXOTH); +#endif + if (res) { + std::cout << "Cannot open output directory \"" << out_dir << "\"" << std::endl; + } + return; + } + closedir(out_dir_p); +} + +void find_dir_and_file_name(std::string a, std::string & dir, std::string& fn) { + // todo: make it system independent + size_t last_slash_pos = a.find_last_of("/"); + if (last_slash_pos >= a.size()) { + std::cout << "cannot find file name in " << a << std::endl; + throw; + } + dir = a.substr(0, last_slash_pos); + // std::cout << "dir = " << dir << std::endl; + fn = a.substr(last_slash_pos + 1); + // std::cout << "fn = " << fn << std::endl; +} + +void process_test_file(std::string test_dir, std::string test_file_name, argument_parser & args_parser, std::string out_dir, unsigned max_iters, unsigned time_limit, unsigned & successes, unsigned & failures, unsigned & inconclusives); + +void solve_some_mps(argument_parser & args_parser) { + unsigned max_iters, time_limit; + get_time_limit_and_max_iters_from_parser(args_parser, time_limit, max_iters); + unsigned successes = 0; + unsigned failures = 0; + unsigned inconclusives = 0; + std::set minimums; + vector file_names; + fill_file_names(file_names, minimums); + bool solve_for_rational = args_parser.option_is_used("--mpq"); + bool dual = args_parser.option_is_used("--dual"); + bool compare_with_primal = args_parser.option_is_used("--compare_with_primal"); + bool compare_with_glpk = args_parser.option_is_used("--compare_with_glpk"); + if (compare_with_glpk) { + std::string out_dir = args_parser.get_option_value("--out_dir"); + if (out_dir.size() == 0) { + out_dir = "/tmp/test"; + } + test_out_dir(out_dir); + for (auto& a : file_names) { + try { + std::string file_dir; + std::string file_name; + find_dir_and_file_name(a, file_dir, file_name); + process_test_file(file_dir, file_name, args_parser, out_dir, max_iters, time_limit, successes, failures, inconclusives); + } + catch(const char *s){ + std::cout<< "exception: "<< s << std::endl; + } + } + std::cout << "comparing with glpk: successes " << successes << ", failures " << failures << ", inconclusives " << inconclusives << std::endl; + return; + } + if (!solve_for_rational) { + solve_mps(file_names[6], false, 0, time_limit, false, dual, compare_with_primal, args_parser); + solve_mps_with_known_solution(file_names[3], nullptr, INFEASIBLE, dual); // chvatal: 135(d) + std::unordered_map sol; + sol["X1"] = 0; + sol["X2"] = 6; + sol["X3"] = 0; + sol["X4"] = 15; + sol["X5"] = 2; + sol["X6"] = 1; + sol["X7"] = 1; + sol["X8"] = 0; + solve_mps_with_known_solution(file_names[9], &sol, OPTIMAL, dual); + solve_mps_with_known_solution(file_names[0], &sol, OPTIMAL, dual); + sol.clear(); + sol["X1"] = 25.0/14.0; + // sol["X2"] = 0; + // sol["X3"] = 0; + // sol["X4"] = 0; + // sol["X5"] = 0; + // sol["X6"] = 0; + // sol["X7"] = 9.0/14.0; + solve_mps_with_known_solution(file_names[5], &sol, OPTIMAL, dual); // chvatal: 135(e) + solve_mps_with_known_solution(file_names[4], &sol, OPTIMAL, dual); // chvatal: 135(e) + solve_mps_with_known_solution(file_names[2], nullptr, UNBOUNDED, dual); // chvatal: 135(c) + solve_mps_with_known_solution(file_names[1], nullptr, UNBOUNDED, dual); // chvatal: 135(b) + solve_mps(file_names[8], false, 0, time_limit, false, dual, compare_with_primal, args_parser); + // return; + for (auto& s : file_names) { + try { + solve_mps(s, minimums.find(s) != minimums.end(), max_iters, time_limit, false, dual, compare_with_primal, args_parser); + } + catch(const char *s){ + std::cout<< "exception: "<< s << std::endl; + } + } + } else { + // unsigned i = 0; + for (auto& s : file_names) { + // if (i++ > 9) return; + try { + solve_mps_in_rational(s, dual, args_parser); + } + catch(const char *s){ + std::cout<< "exception: "<< s << std::endl; + } + } + } +} +#endif + +void solve_rational() { + lp_primal_simplex solver; + solver.add_constraint(lp_relation::Equal, lean::mpq(7), 0); + solver.add_constraint(lp_relation::Equal, lean::mpq(-3), 1); + + // setting the cost + int cost[] = {-3, -1, -1, 2, -1, 1, 1, -4}; + std::string var_names[8] = {"x1", "x2", "x3", "x4", "x5", "x6", "x7", "x8"}; + + for (unsigned i = 0; i < 8; i++) { + solver.set_cost_for_column(i, lean::mpq(cost[i])); + solver.give_symbolic_name_to_column(var_names[i], i); + } + + int row0[] = {1, 0, 3, 1, -5, -2 , 4, -6}; + for (unsigned i = 0; i < 8; i++) { + solver.set_row_column_coefficient(0, i, lean::mpq(row0[i])); + } + + int row1[] = {0, 1, -2, -1, 4, 1, -3, 5}; + for (unsigned i = 0; i < 8; i++) { + solver.set_row_column_coefficient(1, i, lean::mpq(row1[i])); + } + + int bounds[] = {8, 6, 4, 15, 2, 10, 10, 3}; + for (unsigned i = 0; i < 8; i++) { + solver.set_low_bound(i, lean::mpq(0)); + solver.set_upper_bound(i, lean::mpq(bounds[i])); + } + + std::unordered_map expected_sol; + expected_sol["x1"] = lean::mpq(0); + expected_sol["x2"] = lean::mpq(6); + expected_sol["x3"] = lean::mpq(0); + expected_sol["x4"] = lean::mpq(15); + expected_sol["x5"] = lean::mpq(2); + expected_sol["x6"] = lean::mpq(1); + expected_sol["x7"] = lean::mpq(1); + expected_sol["x8"] = lean::mpq(0); + solver.find_maximal_solution(); + lean_assert(solver.get_status() == OPTIMAL); + for (auto it : expected_sol) { + lean_assert(it.second == solver.get_column_value_by_name(it.first)); + } +} + + +std::string read_line(bool & end, std::ifstream & file) { + std::string s; + if (!getline(file, s)) { + end = true; + return std::string(); + } + end = false; + return s; +} + +bool contains(std::string const & s, char const * pattern) { + return s.find(pattern) != std::string::npos; +} + + +std::unordered_map * get_solution_from_glpsol_output(std::string & file_name) { + std::ifstream file(file_name); + if (!file.is_open()){ + std::cerr << "cannot open " << file_name << std::endl; + return nullptr; + } + std::string s; + bool end; + do { + s = read_line(end, file); + if (end) { + std::cerr << "unexpected file end " << file_name << std::endl; + return nullptr; + } + if (contains(s, "Column name")){ + break; + } + } while (true); + + read_line(end, file); + if (end) { + std::cerr << "unexpected file end " << file_name << std::endl; + return nullptr; + } + + auto ret = new std::unordered_map(); + + do { + s = read_line(end, file); + if (end) { + std::cerr << "unexpected file end " << file_name << std::endl; + return nullptr; + } + auto split = string_split(s, " \t", false); + if (split.size() == 0) { + return ret; + } + + lean_assert(split.size() > 3); + (*ret)[split[1]] = atof(split[3].c_str()); + } while (true); +} + + + +void test_init_U() { + static_matrix m(3, 7); + m(0, 0) = 10; m(0, 1) = 11; m(0, 2) = 12; m(0, 3) = 13; m(0, 4) = 14; + m(1, 0) = 20; m(1, 1) = 21; m(1, 2) = 22; m(1, 3) = 23; m(1, 5) = 24; + m(2, 0) = 30; m(2, 1) = 31; m(2, 2) = 32; m(2, 3) = 33; m(2, 6) = 34; +#ifdef LEAN_DEBUG + print_matrix(m, std::cout); +#endif + vector basis(3); + basis[0] = 1; + basis[1] = 2; + basis[2] = 4; + + sparse_matrix u(m, basis); + + for (unsigned i = 0; i < 3; i++) { + for (unsigned j = 0; j < 3; j ++) { + lean_assert(m(i, basis[j]) == u(i, j)); + } + } + + // print_matrix(m); + // print_matrix(u); +} + +void test_replace_column() { + sparse_matrix m(10); + fill_matrix(m); + m.swap_columns(0, 7); + m.swap_columns(6, 3); + m.swap_rows(2, 0); + for (unsigned i = 1; i < m.dimension(); i++) { + m(i, 0) = 0; + } + + indexed_vector w(m.dimension()); + for (unsigned i = 0; i < m.dimension(); i++) { + w.set_value(i % 3, i); + } + + lp_settings settings; + + for (unsigned column_to_replace = 0; column_to_replace < m.dimension(); column_to_replace ++) { + m.replace_column(column_to_replace, w, settings); + for (unsigned i = 0; i < m.dimension(); i++) { + lean_assert(abs(w[i] - m(i, column_to_replace)) < 0.00000001); + } + } +} + + +void setup_args_parser(argument_parser & parser) { + parser.add_option_with_help_string("-xyz_sample", "run a small interactive scenario"); + parser.add_option_with_after_string_with_help("--density", "the percentage of non-zeroes in the matrix below which it is not dense"); + parser.add_option_with_after_string_with_help("--harris_toler", "harris tolerance"); + parser.add_option_with_help_string("--test_swaps", "test row swaps with a permutation"); + parser.add_option_with_help_string("--test_perm", "test permutaions"); + parser.add_option_with_after_string_with_help("--checklu", "the file name for lu checking"); + parser.add_option_with_after_string_with_help("--partial_pivot", "the partial pivot constant, a number somewhere between 10 and 100"); + parser.add_option_with_after_string_with_help("--percent_for_enter", "which percent of columns check for entering column"); + parser.add_option_with_help_string("--totalinf", "minimizes the total infeasibility instead of diminishin infeasibility of the rows"); + parser.add_option_with_after_string_with_help("--rep_frq", "the report frequency, in how many iterations print the cost and other info "); + parser.add_option_with_help_string("--smt", "smt file format"); + parser.add_option_with_after_string_with_help("--filelist", "the file containing the list of files"); + parser.add_option_with_after_string_with_help("--file", "the input file name"); + parser.add_option_with_after_string_with_help("--random_seed", "random seed"); + parser.add_option_with_help_string("--bp", "bound propogation"); + parser.add_option_with_help_string("--min", "will look for the minimum for the given file if --file is used; the default is looking for the max"); + parser.add_option_with_help_string("--max", "will look for the maximum for the given file if --file is used; it is the default behavior"); + parser.add_option_with_after_string_with_help("--max_iters", "maximum total iterations in a core solver stage"); + parser.add_option_with_after_string_with_help("--time_limit", "time limit in seconds"); + parser.add_option_with_help_string("--mpq", "solve for rational numbers"); + parser.add_option_with_after_string_with_help("--simplex_strategy", "sets simplex strategy for rational number"); + parser.add_option_with_help_string("--test_lu", "test the work of the factorization"); + parser.add_option_with_help_string("--test_small_lu", "test the work of the factorization on a smallish matrix"); + parser.add_option_with_help_string("--test_larger_lu", "test the work of the factorization"); + parser.add_option_with_help_string("--test_larger_lu_with_holes", "test the work of the factorization"); + parser.add_option_with_help_string("--test_lp_0", "solve a small lp"); + parser.add_option_with_help_string("--solve_some_mps", "solves a list of mps problems"); + parser.add_option_with_after_string_with_help("--test_file_directory", "loads files from the directory for testing"); + parser.add_option_with_help_string("--compare_with_glpk", "compares the results by running glpsol"); + parser.add_option_with_after_string_with_help("--out_dir", "setting the output directory for tests, if not set /tmp is used"); + parser.add_option_with_help_string("--dual", "using the dual simplex solver"); + parser.add_option_with_help_string("--compare_with_primal", "using the primal simplex solver for comparison"); + parser.add_option_with_help_string("--lar", "test lar_solver"); + parser.add_option_with_after_string_with_help("--maxng", "max iterations without progress"); + parser.add_option_with_help_string("-tbq", "test binary queue"); + parser.add_option_with_help_string("--randomize_lar", "test randomize funclionality"); + parser.add_option_with_help_string("--smap", "test stacked_map"); + parser.add_option_with_help_string("--term", "simple term test"); + parser.add_option_with_help_string("--eti"," run a small evidence test for total infeasibility scenario"); + parser.add_option_with_help_string("--row_inf", "forces row infeasibility search"); + parser.add_option_with_help_string("-pd", "presolve with double solver"); + parser.add_option_with_help_string("--test_int_set", "test int_set"); + parser.add_option_with_help_string("--test_mpq", "test rationals"); + parser.add_option_with_help_string("--test_mpq_np", "test rationals"); + parser.add_option_with_help_string("--test_mpq_np_plus", "test rationals using plus instead of +="); +} + +struct fff { int a; int b;}; + +void test_stacked_map_itself() { + vector v(3,0); + for(auto u : v) + std::cout << u << std::endl; + + std::unordered_map foo; + fff l; + l.a = 0; + l.b =1; + foo[1] = l; + int r = 1; + int k = foo[r].a; + std::cout << k << std::endl; + + stacked_map m; + m[0] = 3; + m[1] = 4; + m.push(); + m[1] = 5; + m[2] = 2; + m.pop(); + m.erase(2); + m[2] = 3; + m.erase(1); + m.push(); + m[3] = 100; + m[4] = 200; + m.erase(1); + m.push(); + m[5] = 300; + m[6] = 400; + m[5] = 301; + m.erase(5); + m[3] = 122; + + m.pop(2); + m.pop(); +} + +void test_stacked_unsigned() { + std::cout << "test stacked unsigned" << std::endl; + stacked_value v(0); + v = 1; + v = 2; + v.push(); + v = 3; + v = 4; + v.pop(); + lean_assert(v == 2); + v ++; + v++; + std::cout << "before push v=" << v << std::endl; + v.push(); + v++; + v.push(); + v+=1; + std::cout << "v = " << v << std::endl; + v.pop(2); + lean_assert(v == 4); + const unsigned & rr = v; + std::cout << rr << std:: endl; + +} + +void test_stacked_value() { + test_stacked_unsigned(); +} + +void test_stacked_vector() { + std::cout << "test_stacked_vector" << std::endl; + stacked_vector v; + v.push(); + v.push_back(0); + v.push_back(1); + v.push(); + v[0] = 3; + v[0] = 0; + v.push_back(2); + v.push_back(3); + v.push_back(34); + v.push(); + v[1]=3; + v[2] = 3; + v.push(); + v[0]= 7; + v[1] = 9; + v.pop(2); + if (v.size()) + v[v.size() -1 ] = 7; + v.push(); + v.push_back(33); + v[0] = 13; + v.pop(); + +} + +void test_stacked_set() { +#ifdef LEAN_DEBUG + std::cout << "test_stacked_set" << std::endl; + stacked_unordered_set s; + s.insert(1); + s.insert(2); + s.insert(3); + std::unordered_set scopy = s(); + s.push(); + s.insert(4); + s.pop(); + lean_assert(s() == scopy); + s.push(); + s.push(); + s.insert(4); + s.insert(5); + s.push(); + s.insert(4); + s.pop(3); + lean_assert(s() == scopy); +#endif +} + +void test_stacked() { + std::cout << "test_stacked_map()" << std::endl; + test_stacked_map_itself(); + test_stacked_value(); + test_stacked_vector(); + test_stacked_set(); + +} + +char * find_home_dir() { + #ifdef _WINDOWS + #else + char * home_dir = getenv("HOME"); + if (home_dir == nullptr) { + std::cout << "cannot find home directory" << std::endl; + return nullptr; + } + #endif + return nullptr; +} + + +template +void print_chunk(T * arr, unsigned len) { + for (unsigned i = 0; i < len; i++) { + std::cout << arr[i] << ", "; + } + std::cout << std::endl; +} + +struct mem_cpy_place_holder { + static void mem_copy_hook(int * destination, unsigned num) { + if (destination == nullptr || num == 0) { + throw "bad parameters"; + } + } +}; + +void finalize(unsigned ret) { + /* + finalize_util_module(); + finalize_numerics_module(); + */ + // return ret; +} + +void get_time_limit_and_max_iters_from_parser(argument_parser & args_parser, unsigned & time_limit, unsigned & max_iters) { + std::string s = args_parser.get_option_value("--max_iters"); + if (s.size() > 0) { + max_iters = atoi(s.c_str()); + } else { + max_iters = 0; + } + + std::string time_limit_string = args_parser.get_option_value("--time_limit"); + if (time_limit_string.size() > 0) { + time_limit = atoi(time_limit_string.c_str()); + } else { + time_limit = 0; + } +} + + +std::string create_output_file_name(bool minimize, std::string file_name, bool use_mpq) { + std::string ret = file_name + "_lp_tst_" + (minimize?"min":"max"); + if (use_mpq) return ret + "_mpq.out"; + return ret + ".out"; +} + +std::string create_output_file_name_for_glpsol(bool minimize, std::string file_name){ + return file_name + (minimize?"_min":"_max") + "_glpk_out"; +} + +int run_glpk(std::string file_name, std::string glpk_out_file_name, bool minimize, unsigned time_limit) { + std::string minmax(minimize?"--min":"--max"); + std::string tmlim = time_limit > 0 ? std::string(" --tmlim ") + std::to_string(time_limit)+ " ":std::string(); + std::string command_line = std::string("glpsol --nointopt --nomip ") + minmax + tmlim + + " -o " + glpk_out_file_name +" " + file_name + " > /dev/null"; + return system(command_line.c_str()); +} + +std::string get_status(std::string file_name) { + std::ifstream f(file_name); + if (!f.is_open()) { + std::cout << "cannot open " << file_name << std::endl; + throw 0; + } + std::string str; + while (getline(f, str)) { + if (str.find("Status") != std::string::npos) { + vector tokens = split_and_trim(str); + if (tokens.size() != 2) { + std::cout << "unexpected Status string " << str << std::endl; + throw 0; + } + return tokens[1]; + } + } + std::cout << "cannot find the status line in " << file_name << std::endl; + throw 0; +} + +// returns true if the costs should be compared too +bool compare_statuses(std::string glpk_out_file_name, std::string lp_out_file_name, unsigned & successes, unsigned & failures) { + std::string glpk_status = get_status(glpk_out_file_name); + std::string lp_tst_status = get_status(lp_out_file_name); + + if (glpk_status != lp_tst_status) { + if (glpk_status == "UNDEFINED" && (lp_tst_status == "UNBOUNDED" || lp_tst_status == "INFEASIBLE")) { + successes++; + return false; + } else { + std::cout << "glpsol and lp_tst disagree: glpsol status is " << glpk_status; + std::cout << " but lp_tst status is " << lp_tst_status << std::endl; + failures++; + return false; + } + } + return lp_tst_status == "OPTIMAL"; +} + +double get_glpk_cost(std::string file_name) { + std::ifstream f(file_name); + if (!f.is_open()) { + std::cout << "cannot open " << file_name << std::endl; + throw 0; + } + std::string str; + while (getline(f, str)) { + if (str.find("Objective") != std::string::npos) { + vector tokens = split_and_trim(str); + if (tokens.size() != 5) { + std::cout << "unexpected Objective std::string " << str << std::endl; + throw 0; + } + return atof(tokens[3].c_str()); + } + } + std::cout << "cannot find the Objective line in " << file_name << std::endl; + throw 0; +} + +double get_lp_tst_cost(std::string file_name) { + std::ifstream f(file_name); + if (!f.is_open()) { + std::cout << "cannot open " << file_name << std::endl; + throw 0; + } + std::string str; + std::string cost_string; + while (getline(f, str)) { + if (str.find("cost") != std::string::npos) { + cost_string = str; + } + } + if (cost_string.size() == 0) { + std::cout << "cannot find the cost line in " << file_name << std::endl; + throw 0; + } + + vector tokens = split_and_trim(cost_string); + if (tokens.size() != 3) { + std::cout << "unexpected cost string " << cost_string << std::endl; + throw 0; + } + return atof(tokens[2].c_str()); +} + +bool values_are_one_percent_close(double a, double b) { + double maxval = std::max(fabs(a), fabs(b)); + if (maxval < 0.000001) { + return true; + } + + double one_percent = maxval / 100; + return fabs(a - b) <= one_percent; +} + +// returns true if both are optimal +void compare_costs(std::string glpk_out_file_name, + std::string lp_out_file_name, + unsigned & successes, + unsigned & failures) { + double a = get_glpk_cost(glpk_out_file_name); + double b = get_lp_tst_cost(lp_out_file_name); + + if (values_are_one_percent_close(a, b)) { + successes++; + } else { + failures++; + std::cout << "glpsol cost is " << a << " lp_tst cost is " << b << std::endl; + } +} + + + +void compare_with_glpk(std::string glpk_out_file_name, std::string lp_out_file_name, unsigned & successes, unsigned & failures, std::string /*lp_file_name*/) { +#ifdef CHECK_GLPK_SOLUTION + std::unordered_map * solution_table = get_solution_from_glpsol_output(glpk_out_file_name); + if (solution_is_feasible(lp_file_name, *solution_table)) { + std::cout << "glpk solution is feasible" << std::endl; + } else { + std::cout << "glpk solution is infeasible" << std::endl; + } + delete solution_table; +#endif + if (compare_statuses(glpk_out_file_name, lp_out_file_name, successes, failures)) { + compare_costs(glpk_out_file_name, lp_out_file_name, successes, failures); + } +} +void test_lar_on_file(std::string file_name, argument_parser & args_parser); + +void process_test_file(std::string test_dir, std::string test_file_name, argument_parser & args_parser, std::string out_dir, unsigned max_iters, unsigned time_limit, unsigned & successes, unsigned & failures, unsigned & inconclusives) { + bool use_mpq = args_parser.option_is_used("--mpq"); + bool minimize = args_parser.option_is_used("--min"); + std::string full_lp_tst_out_name = out_dir + "/" + create_output_file_name(minimize, test_file_name, use_mpq); + + std::string input_file_name = test_dir + "/" + test_file_name; + if (input_file_name[input_file_name.size() - 1] == '~') { + // std::cout << "ignoring " << input_file_name << std::endl; + return; + } + std::cout <<"processing " << input_file_name << std::endl; + + std::ofstream out(full_lp_tst_out_name); + if (!out.is_open()) { + std::cout << "cannot open file " << full_lp_tst_out_name << std::endl; + throw 0; + } + std::streambuf *coutbuf = std::cout.rdbuf(); // save old buffer + std::cout.rdbuf(out.rdbuf()); // redirect std::cout to dir_entry->d_name! + bool dual = args_parser.option_is_used("--dual"); + try { + if (args_parser.option_is_used("--lar")) + test_lar_on_file(input_file_name, args_parser); + else + solve_mps(input_file_name, minimize, max_iters, time_limit, use_mpq, dual, false, args_parser); + } + catch(...) { + std::cout << "catching the failure" << std::endl; + failures++; + std::cout.rdbuf(coutbuf); // reset to standard output again + return; + } + std::cout.rdbuf(coutbuf); // reset to standard output again + + if (args_parser.option_is_used("--compare_with_glpk")) { + std::string glpk_out_file_name = out_dir + "/" + create_output_file_name_for_glpsol(minimize, std::string(test_file_name)); + int glpk_exit_code = run_glpk(input_file_name, glpk_out_file_name, minimize, time_limit); + if (glpk_exit_code != 0) { + std::cout << "glpk failed" << std::endl; + inconclusives++; + } else { + compare_with_glpk(glpk_out_file_name, full_lp_tst_out_name, successes, failures, input_file_name); + } + } +} +/* + int my_readdir(DIR *dirp, struct dirent * +#ifndef LEAN_WINDOWS + entry +#endif + , struct dirent **result) { +#ifdef LEAN_WINDOWS + *result = readdir(dirp); // NOLINT + return *result != nullptr? 0 : 1; +#else + return readdir_r(dirp, entry, result); +#endif +} +*/ +/* +vector> get_file_list_of_dir(std::string test_file_dir) { + DIR *dir; + if ((dir = opendir(test_file_dir.c_str())) == nullptr) { + std::cout << "Cannot open directory " << test_file_dir << std::endl; + throw 0; + } + vector> ret; + struct dirent entry; + struct dirent* result; + int return_code; + for (return_code = my_readdir(dir, &entry, &result); +#ifndef LEAN_WINDOWS + result != nullptr && +#endif + return_code == 0; + return_code = my_readdir(dir, &entry, &result)) { + DIR *tmp_dp = opendir(result->d_name); + struct stat file_record; + if (tmp_dp == nullptr) { + std::string s = test_file_dir+ "/" + result->d_name; + int stat_ret = stat(s.c_str(), & file_record); + if (stat_ret!= -1) { + ret.push_back(make_pair(result->d_name, file_record.st_size)); + } else { + perror("stat"); + exit(1); + } + } else { + closedir(tmp_dp); + } + } + closedir(dir); + return ret; +} +*/ +/* +struct file_size_comp { + unordered_map& m_file_sizes; + file_size_comp(unordered_map& fs) :m_file_sizes(fs) {} + int operator()(std::string a, std::string b) { + std::cout << m_file_sizes.size() << std::endl; + std::cout << a << std::endl; + std::cout << b << std::endl; + + auto ls = m_file_sizes.find(a); + std::cout << "fa" << std::endl; + auto rs = m_file_sizes.find(b); + std::cout << "fb" << std::endl; + if (ls != m_file_sizes.end() && rs != m_file_sizes.end()) { + std::cout << "fc " << std::endl; + int r = (*ls < *rs? -1: (*ls > *rs)? 1 : 0); + std::cout << "calc r " << std::endl; + return r; + } else { + std::cout << "sc " << std::endl; + return 0; + } + } +}; + +*/ +struct sort_pred { + bool operator()(const std::pair &left, const std::pair &right) { + return left.second < right.second; + } +}; + + +void test_files_from_directory(std::string test_file_dir, argument_parser & args_parser) { + /* + std::cout << "loading files from directory \"" << test_file_dir << "\"" << std::endl; + std::string out_dir = args_parser.get_option_value("--out_dir"); + if (out_dir.size() == 0) { + out_dir = "/tmp/test"; + } + DIR *out_dir_p = opendir(out_dir.c_str()); + if (out_dir_p == nullptr) { + std::cout << "Cannot open output directory \"" << out_dir << "\"" << std::endl; + return; + } + closedir(out_dir_p); + vector> files = get_file_list_of_dir(test_file_dir); + std::sort(files.begin(), files.end(), sort_pred()); + unsigned max_iters, time_limit; + get_time_limit_and_max_iters_from_parser(args_parser, time_limit, max_iters); + unsigned successes = 0, failures = 0, inconclusives = 0; + for (auto & t : files) { + process_test_file(test_file_dir, t.first, args_parser, out_dir, max_iters, time_limit, successes, failures, inconclusives); + } + std::cout << "comparing with glpk: successes " << successes << ", failures " << failures << ", inconclusives " << inconclusives << std::endl; + */ +} + + +std::unordered_map get_solution_map(lp_solver * lps, mps_reader & reader) { + std::unordered_map ret; + for (auto it : reader.column_names()) { + ret[it] = lps->get_column_value_by_name(it); + } + return ret; +} + +void run_lar_solver(argument_parser & args_parser, lar_solver * solver, mps_reader * reader) { + std::string maxng = args_parser.get_option_value("--maxng"); + if (maxng.size() > 0) { + solver->settings().max_number_of_iterations_with_no_improvements = atoi(maxng.c_str()); + } + if (args_parser.option_is_used("-pd")){ + solver->settings().presolve_with_double_solver_for_lar = true; + } + + std::string iter = args_parser.get_option_value("--max_iters"); + if (iter.size() > 0) { + solver->settings().max_total_number_of_iterations = atoi(iter.c_str()); + } + if (args_parser.option_is_used("--compare_with_primal")){ + if (reader == nullptr) { + std::cout << "cannot compare with primal, the reader is null " << std::endl; + return; + } + auto * lps = reader->create_solver(false); + lps->find_maximal_solution(); + std::unordered_map sol = get_solution_map(lps, *reader); + std::cout << "status = " << lp_status_to_string(solver->get_status()) << std::endl; + return; + } + int begin = get_millisecond_count(); + lp_status status = solver->solve(); + std::cout << "status is " << lp_status_to_string(status) << ", processed for " << get_millisecond_span(begin) / 1000.0 <<" seconds, and " << solver->get_total_iterations() << " iterations" << std::endl; + if (solver->get_status() == INFEASIBLE) { + vector> evidence; + solver->get_infeasibility_explanation(evidence); + } + if (args_parser.option_is_used("--randomize_lar")) { + if (solver->get_status() != OPTIMAL) { + std::cout << "cannot check randomize on an infeazible problem" << std::endl; + return; + } + std::cout << "checking randomize" << std::endl; + vector all_vars = solver->get_list_of_all_var_indices(); + unsigned m = all_vars.size(); + if (m > 100) + m = 100; + + var_index *vars = new var_index[m]; + for (unsigned i = 0; i < m; i++) + vars[i]=all_vars[i]; + + solver->random_update(m, vars); + delete []vars; + } +} + +lar_solver * create_lar_solver_from_file(std::string file_name, argument_parser & args_parser) { + if (args_parser.option_is_used("--smt")) { + smt_reader reader(file_name); + reader.read(); + if (!reader.is_ok()){ + std::cout << "cannot process " << file_name << std::endl; + return nullptr; + } + return reader.create_lar_solver(); + } + mps_reader reader(file_name); + reader.read(); + if (!reader.is_ok()) { + std::cout << "cannot process " << file_name << std::endl; + return nullptr; + } + return reader.create_lar_solver(); +} + +void test_lar_on_file(std::string file_name, argument_parser & args_parser) { + lar_solver * solver = create_lar_solver_from_file(file_name, args_parser); + mps_reader reader(file_name); + mps_reader * mps_reader = nullptr; + reader.read(); + if (reader.is_ok()) { + mps_reader = & reader; + run_lar_solver(args_parser, solver, mps_reader); + } + delete solver; +} + +vector get_file_names_from_file_list(std::string filelist) { + std::ifstream file(filelist); + if (!file.is_open()) { + std::cout << "cannot open " << filelist << std::endl; + return vector(); + } + vector ret; + bool end; + do { + std::string s = read_line(end, file); + if (end) + break; + if (s.size() == 0) + break; + ret.push_back(s); + } while (true); + return ret; +} + +void test_lar_solver(argument_parser & args_parser) { + + std::string file_name = args_parser.get_option_value("--file"); + if (file_name.size() > 0) { + test_lar_on_file(file_name, args_parser); + return; + } + + std::string file_list = args_parser.get_option_value("--filelist"); + if (file_list.size() > 0) { + for (std::string fn : get_file_names_from_file_list(file_list)) + test_lar_on_file(fn, args_parser); + return; + } +} + +void test_numeric_pair() { + numeric_pair a; + numeric_pair b(2, lean::mpq(6, 2)); + a = b; + numeric_pair c(0.1, 0.5); + a += 2*c; + a -= c; + lean_assert (a == b + c); + numeric_pair d = a * 2; + std::cout << a << std::endl; + lean_assert(b == b); + lean_assert(b < a); + lean_assert(b <= a); + lean_assert(a > b); + lean_assert(a != b); + lean_assert(a >= b); + lean_assert(-a < b); + lean_assert(a < 2 * b); + lean_assert(b + b > a); + lean_assert(lean::mpq(2.1) * b + b > a); + lean_assert(-b * lean::mpq(2.1) - b < lean::mpq(0.99) * a); + std::cout << - b * lean::mpq(2.1) - b << std::endl; + lean_assert(-b *(lean::mpq(2.1) + 1) == - b * lean::mpq(2.1) - b); +} + +void get_matrix_dimensions(std::ifstream & f, unsigned & m, unsigned & n) { + std::string line; + getline(f, line); + getline(f, line); + vector r = split_and_trim(line); + m = atoi(r[1].c_str()); + getline(f, line); + r = split_and_trim(line); + n = atoi(r[1].c_str()); +} + +void read_row_cols(unsigned i, static_matrix& A, std::ifstream & f) { + do { + std::string line; + getline(f, line); + if (line== "row_end") + break; + auto r = split_and_trim(line); + lean_assert(r.size() == 4); + unsigned j = atoi(r[1].c_str()); + double v = atof(r[3].c_str()); + A.set(i, j, v); + } while (true); +} + +bool read_row(static_matrix & A, std::ifstream & f) { + std::string line; + getline(f, line); + if (static_cast(line.find("row")) == -1) + return false; + auto r = split_and_trim(line); + if (r[0] != "row") + std::cout << "wrong row line" << line << std::endl; + unsigned i = atoi(r[1].c_str()); + read_row_cols(i, A, f); + return true; +} + +void read_rows(static_matrix& A, std::ifstream & f) { + while (read_row(A, f)) {} +} + +void read_basis(vector & basis, std::ifstream & f) { + std::cout << "reading basis" << std::endl; + std::string line; + getline(f, line); + lean_assert(line == "basis_start"); + do { + getline(f, line); + if (line == "basis_end") + break; + unsigned j = atoi(line.c_str()); + basis.push_back(j); + } while (true); +} + +void read_indexed_vector(indexed_vector & v, std::ifstream & f) { + std::string line; + getline(f, line); + lean_assert(line == "vector_start"); + do { + getline(f, line); + if (line == "vector_end") break; + auto r = split_and_trim(line); + unsigned i = atoi(r[0].c_str()); + double val = atof(r[1].c_str()); + v.set_value(val, i); + std::cout << "setting value " << i << " = " << val << std::endl; + } while (true); +} + +void check_lu_from_file(std::string lufile_name) { + std::ifstream f(lufile_name); + if (!f.is_open()) { + std::cout << "cannot open file " << lufile_name << std::endl; + } + unsigned m, n; + get_matrix_dimensions(f, m, n); + std::cout << "init matrix " << m << " by " << n << std::endl; + static_matrix A(m, n); + read_rows(A, f); + vector basis; + read_basis(basis, f); + indexed_vector v(m); + // read_indexed_vector(v, f); + f.close(); + vector basis_heading; + lp_settings settings; + vector non_basic_columns; + lu lsuhl(A, basis, settings); + indexed_vector d(A.row_count()); + unsigned entering = 26; + lsuhl.solve_Bd(entering, d, v); +#ifdef LEAN_DEBUG + auto B = get_B(lsuhl, basis); + vector a(m); + A.copy_column_to_vector(entering, a); + indexed_vector cd(d); + B.apply_from_left(cd.m_data, settings); + lean_assert(vectors_are_equal(cd.m_data , a)); +#endif +} + +void test_square_dense_submatrix() { + std::cout << "testing square_dense_submatrix" << std::endl; + unsigned parent_dim = 7; + sparse_matrix parent(parent_dim); + fill_matrix(parent); + unsigned index_start = 3; + square_dense_submatrix d; + d.init(&parent, index_start); + for (unsigned i = index_start; i < parent_dim; i++) + for (unsigned j = index_start; j < parent_dim; j++) + d[i][j] = i*3+j*2; +#ifdef LEAN_DEBUG + unsigned dim = parent_dim - index_start; + dense_matrix m(dim, dim); + for (unsigned i = index_start; i < parent_dim; i++) + for (unsigned j = index_start; j < parent_dim; j++) + m[i-index_start][j-index_start] = d[i][j]; + print_matrix(&m, std::cout); +#endif + for (unsigned i = index_start; i < parent_dim; i++) + for (unsigned j = index_start; j < parent_dim; j++) + d[i][j] = d[j][i]; +#ifdef LEAN_DEBUG + for (unsigned i = index_start; i < parent_dim; i++) + for (unsigned j = index_start; j < parent_dim; j++) + m[i-index_start][j-index_start] = d[i][j]; + + print_matrix(&m, std::cout); + std::cout << std::endl; +#endif +} + + + +void print_st(lp_status status) { + std::cout << lp_status_to_string(status) << std::endl; +} + + + +void test_term() { + lar_solver solver; + unsigned _x = 0; + unsigned _y = 1; + var_index x = solver.add_var(_x); + var_index y = solver.add_var(_y); + + vector> term_ls; + term_ls.push_back(std::pair((int)1, x)); + term_ls.push_back(std::pair((int)1, y)); + var_index z = solver.add_term(term_ls, mpq(3)); + + vector> ls; + ls.push_back(std::pair((int)1, x)); + ls.push_back(std::pair((int)1, y)); + ls.push_back(std::pair((int)1, z)); + + solver.add_constraint(ls, lconstraint_kind::EQ, mpq(0)); + auto status = solver.solve(); + std::cout << lp_status_to_string(status) << std::endl; + std::unordered_map model; + solver.get_model(model); + + for (auto & t : model) { + std::cout << solver.get_variable_name(t.first) << " = " << t.second.get_double() << ","; + } + std::cout << std::endl; + +} + +void test_evidence_for_total_inf_simple(argument_parser & args_parser) { + lar_solver solver; + var_index x = solver.add_var(0); + var_index y = solver.add_var(1); + solver.add_var_bound(x, LE, -mpq(1)); + solver.add_var_bound(y, GE, mpq(0)); + vector> ls; + + ls.push_back(std::pair((int)1, x)); + ls.push_back(std::pair((int)1, y)); + solver.add_constraint(ls, GE, mpq(1)); + ls.pop_back(); + ls.push_back(std::pair(-(int)1, y)); + solver.add_constraint(ls, lconstraint_kind::GE, mpq(0)); + auto status = solver.solve(); + std::cout << lp_status_to_string(status) << std::endl; + std::unordered_map model; + lean_assert(solver.get_status() == INFEASIBLE); +} +void test_bound_propagation_one_small_sample1() { + /* +(<= (+ a (* (- 1.0) b)) 0.0) +(<= (+ b (* (- 1.0) x_13)) 0.0) +--> (<= (+ a (* (- 1.0) c)) 0.0) + +the inequality on (<= a c) is obtained from a triangle inequality (<= a b) (<= b c). +If b becomes basic variable, then it is likely the old solver ends up with a row that implies (<= a c). + a - b <= 0.0 + b - c <= 0.0 + + got to get a <= c + */ + std::function bound_is_relevant = + [&](unsigned j, bool is_low_bound, bool strict, const rational& bound_val) { + return true; + }; + lar_solver ls; + unsigned a = ls.add_var(0); + unsigned b = ls.add_var(1); + unsigned c = ls.add_var(2); + vector> coeffs; + coeffs.push_back(std::pair(1, a)); + coeffs.push_back(std::pair(-1, c)); + ls.add_term(coeffs, zero_of_type()); + coeffs.pop_back(); + coeffs.push_back(std::pair(-1, b)); + ls.add_term(coeffs, zero_of_type()); + coeffs.clear(); + coeffs.push_back(std::pair(1, a)); + coeffs.push_back(std::pair(-1, b)); + ls.add_constraint(coeffs, LE, zero_of_type()); + coeffs.clear(); + coeffs.push_back(std::pair(1, b)); + coeffs.push_back(std::pair(-1, c)); + ls.add_constraint(coeffs, LE, zero_of_type()); + vector ev; + ls.add_var_bound(a, LE, mpq(1)); + ls.solve(); + bound_propagator bp(ls); + ls.propagate_bounds_for_touched_rows(bp); + std::cout << " bound ev from test_bound_propagation_one_small_sample1" << std::endl; + for (auto & be : bp.m_ibounds) { + std::cout << "bound\n"; + ls.print_implied_bound(be, std::cout); + } +} + +void test_bound_propagation_one_small_samples() { + test_bound_propagation_one_small_sample1(); + /* + (>= x_46 0.0) +(<= x_29 0.0) +(not (<= x_68 0.0)) +(<= (+ (* (/ 1001.0 1998.0) x_10) (* (- 1.0) x_151) x_68) (- (/ 1001.0 999.0))) +(<= (+ (* (/ 1001.0 999.0) x_9) + (* (- 1.0) x_152) + (* (/ 1001.0 999.0) x_151) + (* (/ 1001.0 999.0) x_68)) + (- (/ 1502501.0 999000.0))) +(not (<= (+ (* (/ 999.0 2.0) x_10) (* (- 1.0) x_152) (* (- (/ 999.0 2.0)) x_151)) + (/ 1001.0 2.0))) +(not (<= x_153 0.0))z +(>= (+ x_9 (* (- (/ 1001.0 999.0)) x_10) (* (- 1.0) x_153) (* (- 1.0) x_68)) + (/ 5003.0 1998.0)) +--> (not (<= (+ x_10 x_46 (* (- 1.0) x_29)) 0.0)) + +and + +(<= (+ a (* (- 1.0) b)) 0.0) +(<= (+ b (* (- 1.0) x_13)) 0.0) +--> (<= (+ a (* (- 1.0) x_13)) 0.0) + +In the first case, there typically are no atomic formulas for bounding x_10. So there is never some +basic lemma of the form (>= x46 0), (<= x29 0), (>= x10 0) -> (not (<= (+ x10 x46 (- x29)) 0)). +Instead the bound on x_10 falls out from a bigger blob of constraints. + +In the second case, the inequality on (<= x19 x13) is obtained from a triangle inequality (<= x19 x9) (<= x9 x13). +If x9 becomes basic variable, then it is likely the old solver ends up with a row that implies (<= x19 x13). + */ +} +void test_bound_propagation_one_row() { + lar_solver ls; + unsigned x0 = ls.add_var(0); + unsigned x1 = ls.add_var(1); + vector> c; + c.push_back(std::pair(1, x0)); + c.push_back(std::pair(-1, x1)); + ls.add_constraint(c, EQ, one_of_type()); + vector ev; + ls.add_var_bound(x0, LE, mpq(1)); + ls.solve(); + bound_propagator bp(ls); + ls.propagate_bounds_for_touched_rows(bp); +} +void test_bound_propagation_one_row_with_bounded_vars() { + lar_solver ls; + unsigned x0 = ls.add_var(0); + unsigned x1 = ls.add_var(1); + vector> c; + c.push_back(std::pair(1, x0)); + c.push_back(std::pair(-1, x1)); + ls.add_constraint(c, EQ, one_of_type()); + vector ev; + ls.add_var_bound(x0, GE, mpq(-3)); + ls.add_var_bound(x0, LE, mpq(3)); + ls.add_var_bound(x0, LE, mpq(1)); + ls.solve(); + bound_propagator bp(ls); + ls.propagate_bounds_for_touched_rows(bp); +} +void test_bound_propagation_one_row_mixed() { + lar_solver ls; + unsigned x0 = ls.add_var(0); + unsigned x1 = ls.add_var(1); + vector> c; + c.push_back(std::pair(1, x0)); + c.push_back(std::pair(-1, x1)); + ls.add_constraint(c, EQ, one_of_type()); + vector ev; + ls.add_var_bound(x1, LE, mpq(1)); + ls.solve(); + bound_propagator bp(ls); + ls.propagate_bounds_for_touched_rows(bp); +} + +void test_bound_propagation_two_rows() { + lar_solver ls; + unsigned x = ls.add_var(0); + unsigned y = ls.add_var(1); + unsigned z = ls.add_var(2); + vector> c; + c.push_back(std::pair(1, x)); + c.push_back(std::pair(2, y)); + c.push_back(std::pair(3, z)); + ls.add_constraint(c, GE, one_of_type()); + c.clear(); + c.push_back(std::pair(3, x)); + c.push_back(std::pair(2, y)); + c.push_back(std::pair(1, z)); + ls.add_constraint(c, GE, one_of_type()); + ls.add_var_bound(x, LE, mpq(2)); + vector ev; + ls.add_var_bound(y, LE, mpq(1)); + ls.solve(); + bound_propagator bp(ls); + ls.propagate_bounds_for_touched_rows(bp); +} + +void test_total_case_u() { + std::cout << "test_total_case_u\n"; + lar_solver ls; + unsigned x = ls.add_var(0); + unsigned y = ls.add_var(1); + unsigned z = ls.add_var(2); + vector> c; + c.push_back(std::pair(1, x)); + c.push_back(std::pair(2, y)); + c.push_back(std::pair(3, z)); + ls.add_constraint(c, LE, one_of_type()); + ls.add_var_bound(x, GE, zero_of_type()); + ls.add_var_bound(y, GE, zero_of_type()); + vector ev; + ls.add_var_bound(z, GE, zero_of_type()); + ls.solve(); + bound_propagator bp(ls); + ls.propagate_bounds_for_touched_rows(bp); +} +bool contains_j_kind(unsigned j, lconstraint_kind kind, const mpq & rs, const vector & ev) { + for (auto & e : ev) { + if (e.m_j == j && e.m_bound == rs && e.kind() == kind) + return true; + } + return false; +} +void test_total_case_l(){ + std::cout << "test_total_case_l\n"; + lar_solver ls; + unsigned x = ls.add_var(0); + unsigned y = ls.add_var(1); + unsigned z = ls.add_var(2); + vector> c; + c.push_back(std::pair(1, x)); + c.push_back(std::pair(2, y)); + c.push_back(std::pair(3, z)); + ls.add_constraint(c, GE, one_of_type()); + ls.add_var_bound(x, LE, one_of_type()); + ls.add_var_bound(y, LE, one_of_type()); + ls.settings().presolve_with_double_solver_for_lar = true; + vector ev; + ls.add_var_bound(z, LE, zero_of_type()); + ls.solve(); + bound_propagator bp(ls); + ls.propagate_bounds_for_touched_rows(bp); + lean_assert(ev.size() == 4); + lean_assert(contains_j_kind(x, GE, - one_of_type(), ev)); +} +void test_bound_propagation() { + test_total_case_u(); + test_bound_propagation_one_small_samples(); + test_bound_propagation_one_row(); + test_bound_propagation_one_row_with_bounded_vars(); + test_bound_propagation_two_rows(); + test_bound_propagation_one_row_mixed(); + test_total_case_l(); + +} + +void test_int_set() { + int_set s(4); + s.insert(2); + s.print(std::cout); + s.insert(1); + s.insert(2); + s.print(std::cout); + lean_assert(s.contains(2)); + lean_assert(s.size() == 2); + s.erase(2); + lean_assert(s.size() == 1); + s.erase(2); + lean_assert(s.size() == 1); + s.print(std::cout); + s.insert(3); + s.insert(2); + s.clear(); + lean_assert(s.size() == 0); + + +} + +void test_rationals_no_numeric_pairs() { + stopwatch sw; + + vector c; + for (unsigned j = 0; j < 10; j ++) + c.push_back(mpq(my_random()%100, 1 + my_random()%100 )); + + vector x; + for (unsigned j = 0; j < 10; j ++) + x.push_back(mpq(my_random()%100, 1 + my_random()%100 )); + + unsigned k = 500000; + mpq r=zero_of_type(); + sw.start(); + + for (unsigned j = 0; j < k; j++){ + mpq val = zero_of_type(); + for (unsigned j=0;j< c.size(); j++){ + val += c[j]*x[j]; + } + + r += val; + } + + sw.stop(); + std::cout << "operation with rationals no pairs " << sw.get_seconds() << std::endl; + std::cout << T_to_string(r) << std::endl; +} + +void test_rationals_no_numeric_pairs_plus() { + stopwatch sw; + + vector c; + for (unsigned j = 0; j < 10; j ++) + c.push_back(mpq(my_random()%100, 1 + my_random()%100 )); + + vector x; + for (unsigned j = 0; j < 10; j ++) + x.push_back(mpq(my_random()%100, 1 + my_random()%100 )); + + unsigned k = 500000; + mpq r=zero_of_type(); + sw.start(); + + for (unsigned j = 0; j < k; j++){ + mpq val = zero_of_type(); + for (unsigned j=0;j< c.size(); j++){ + val = val + c[j]*x[j]; + } + + r = r + val; + } + + sw.stop(); + std::cout << "operation with rationals no pairs " << sw.get_seconds() << std::endl; + std::cout << T_to_string(r) << std::endl; +} + + + +void test_rationals() { + stopwatch sw; + + vector c; + for (unsigned j = 0; j < 10; j ++) + c.push_back(mpq(my_random()%100, 1 + my_random()%100)); + + + + vector> x; + for (unsigned j = 0; j < 10; j ++) + x.push_back(mpq(my_random()%100, 1 + my_random()%100 )); + + std::cout << "x = "; + print_vector(x, std::cout); + + unsigned k = 1000000; + numeric_pair r=zero_of_type>(); + sw.start(); + + for (unsigned j = 0; j < k; j++) { + for (unsigned i = 0; i < c.size(); i++) { + r+= c[i] * x[i]; + } + } + sw.stop(); + std::cout << "operation with rationals " << sw.get_seconds() << std::endl; + std::cout << T_to_string(r) << std::endl; +} + +void test_lp_local(int argn, char**argv) { + std::cout << "resize\n"; + vector r; + r.resize(1); + + // initialize_util_module(); + // initialize_numerics_module(); + int ret; + argument_parser args_parser(argn, argv); + setup_args_parser(args_parser); + if (!args_parser.parse()) { + std::cout << args_parser.m_error_message << std::endl; + std::cout << args_parser.usage_string(); + ret = 1; + return finalize(ret); + } + + args_parser.print(); + + if (args_parser.option_is_used("--test_mpq")) { + test_rationals(); + return finalize(0); + } + + if (args_parser.option_is_used("--test_mpq_np")) { + test_rationals_no_numeric_pairs(); + return finalize(0); + } + + if (args_parser.option_is_used("--test_mpq_np_plus")) { + test_rationals_no_numeric_pairs_plus(); + return finalize(0); + } + + + + if (args_parser.option_is_used("--test_int_set")) { + test_int_set(); + return finalize(0); + } + if (args_parser.option_is_used("--bp")) { + test_bound_propagation(); + return finalize(0); + } + + + std::string lufile = args_parser.get_option_value("--checklu"); + if (lufile.size()) { + check_lu_from_file(lufile); + return finalize(0); + } + +#ifdef LEAN_DEBUG + if (args_parser.option_is_used("--test_swaps")) { + sparse_matrix m(10); + fill_matrix(m); + test_swap_rows_with_permutation(m); + test_swap_cols_with_permutation(m); + return finalize(0); + } +#endif + if (args_parser.option_is_used("--test_perm")) { + test_permutations(); + return finalize(0); + } + if (args_parser.option_is_used("--test_file_directory")) { + test_files_from_directory(args_parser.get_option_value("--test_file_directory"), args_parser); + return finalize(0); + } + std::string file_list = args_parser.get_option_value("--filelist"); + if (file_list.size() > 0) { + for (std::string fn : get_file_names_from_file_list(file_list)) + solve_mps(fn, args_parser); + return finalize(0); + } + + if (args_parser.option_is_used("-tbq")) { + test_binary_priority_queue(); + ret = 0; + return finalize(ret); + } + +#ifdef LEAN_DEBUG + lp_settings settings; + update_settings(args_parser, settings); + if (args_parser.option_is_used("--test_lu")) { + test_lu(settings); + ret = 0; + return finalize(ret); + } + + if (args_parser.option_is_used("--test_small_lu")) { + test_small_lu(settings); + ret = 0; + return finalize(ret); + } + + if (args_parser.option_is_used("--lar")){ + std::cout <<"calling test_lar_solver" << std::endl; + test_lar_solver(args_parser); + return finalize(0); + } + + + + if (args_parser.option_is_used("--test_larger_lu")) { + test_larger_lu(settings); + ret = 0; + return finalize(ret); + } + + if (args_parser.option_is_used("--test_larger_lu_with_holes")) { + test_larger_lu_with_holes(settings); + ret = 0; + return finalize(ret); + } +#endif + if (args_parser.option_is_used("--eti")) { + test_evidence_for_total_inf_simple(args_parser); + ret = 0; + return finalize(ret); + } + + + if (args_parser.option_is_used("--test_lp_0")) { + test_lp_0(); + ret = 0; + return finalize(ret); + } + + if (args_parser.option_is_used("--smap")) { + test_stacked(); + ret = 0; + return finalize(ret); + } + if (args_parser.option_is_used("--term")) { + test_term(); + ret = 0; + return finalize(ret); + } + unsigned max_iters; + unsigned time_limit; + get_time_limit_and_max_iters_from_parser(args_parser, time_limit, max_iters); + bool dual = args_parser.option_is_used("--dual"); + bool solve_for_rational = args_parser.option_is_used("--mpq"); + std::string file_name = args_parser.get_option_value("--file"); + if (file_name.size() > 0) { + solve_mps(file_name, args_parser.option_is_used("--min"), max_iters, time_limit, solve_for_rational, dual, args_parser.option_is_used("--compare_with_primal"), args_parser); + ret = 0; + return finalize(ret); + } + + if (args_parser.option_is_used("--solve_some_mps")) { +#if _LINUX_ + solve_some_mps(args_parser); +#endif + ret = 0; + return finalize(ret); + } + // lean::ccc = 0; + return finalize(0); + test_init_U(); + test_replace_column(); +#ifdef LEAN_DEBUG + sparse_matrix_with_permutaions_test(); + test_dense_matrix(); + test_swap_operations(); + test_permutations(); + test_pivot_like_swaps_and_pivot(); +#endif + tst1(); + std::cout << "done with LP tests\n"; + return finalize(0); // has_violations() ? 1 : 0); +} +} +void tst_lp(char ** argv, int argc, int& i) { + lean::test_lp_local(argc - 2, argv + 2); +} diff --git a/src/test/lp_main.cpp b/src/test/lp_main.cpp new file mode 100644 index 000000000..a301f38c6 --- /dev/null +++ b/src/test/lp_main.cpp @@ -0,0 +1,14 @@ +void gparams_register_modules(){} +void mem_initialize() {} +void mem_finalize() {} +#include "util/rational.h" +namespace lean { +void test_lp_local(int argc, char**argv); +} +int main(int argn, char**argv){ + rational::initialize(); + lean::test_lp_local(argn, argv); + rational::finalize(); + return 0; +} + diff --git a/src/test/smt_reader.h b/src/test/smt_reader.h new file mode 100644 index 000000000..481985ed4 --- /dev/null +++ b/src/test/smt_reader.h @@ -0,0 +1,392 @@ +/* + Copyright (c) 2013 Microsoft Corporation. All rights reserved. + Released under Apache 2.0 license as described in the file LICENSE. + + Author: Lev Nachmanson +*/ + +#pragma once + +// reads an MPS file reperesenting a Mixed Integer Program +#include +#include +#include +#include "util/lp/lp_primal_simplex.h" +#include "util/lp/lp_dual_simplex.h" +#include "util/lp/lar_solver.h" +#include +#include +#include +#include +#include "util/lp/mps_reader.h" +#include "util/lp/ul_pair.h" +#include "util/lp/lar_constraints.h" +#include +#include +namespace lean { + + template + T from_string(const std::string& str) { + std::istringstream ss(str); + T ret; + ss >> ret; + return ret; + } + + class smt_reader { + public: + struct lisp_elem { + std::string m_head; + std::vector m_elems; + void print() { + if (m_elems.size()) { + std::cout << '('; + std::cout << m_head << ' '; + for (auto & el : m_elems) + el.print(); + + std::cout << ')'; + } else { + std::cout << " " << m_head; + } + } + unsigned size() const { return static_cast(m_elems.size()); } + bool is_simple() const { return size() == 0; } + }; + struct formula_constraint { + lconstraint_kind m_kind; + std::vector> m_coeffs; + mpq m_right_side = numeric_traits::zero(); + void add_pair(mpq c, std::string name) { + m_coeffs.push_back(make_pair(c, name)); + } + }; + + lisp_elem m_formula_lisp_elem; + + std::unordered_map m_name_to_var_index; + std::vector m_constraints; + std::string m_file_name; + std::ifstream m_file_stream; + std::string m_line; + bool m_is_OK = true; + unsigned m_line_number = 0; + smt_reader(std::string file_name): + m_file_name(file_name), m_file_stream(file_name) { + } + + void set_error() { + std::cout << "setting error" << std::endl; + m_is_OK = false; + } + + bool is_ok() { + return m_is_OK; + } + + bool prefix(const char * pr) { + return m_line.find(pr) == 0; + } + + int first_separator() { + unsigned blank_pos = static_cast(m_line.find(' ')); + unsigned br_pos = static_cast(m_line.find('(')); + unsigned reverse_br_pos = static_cast(m_line.find(')')); + return std::min(blank_pos, std::min(br_pos, reverse_br_pos)); + } + + void fill_lisp_elem(lisp_elem & lm) { + if (m_line[0] == '(') + fill_nested_elem(lm); + else + fill_simple_elem(lm); + } + + void fill_simple_elem(lisp_elem & lm) { + int separator = first_separator(); + lean_assert(-1 != separator && separator != 0); + lm.m_head = m_line.substr(0, separator); + m_line = m_line.substr(separator); + } + + void fill_nested_elem(lisp_elem & lm) { + lean_assert(m_line[0] == '('); + m_line = m_line.substr(1); + int separator = first_separator(); + lm.m_head = m_line.substr(0, separator); + m_line = m_line.substr(lm.m_head.size()); + eat_blanks(); + while (m_line.size()) { + if (m_line[0] == '(') { + lisp_elem el; + fill_nested_elem(el); + lm.m_elems.push_back(el); + } else { + if (m_line[0] == ')') { + m_line = m_line.substr(1); + break; + } + lisp_elem el; + fill_simple_elem(el); + lm.m_elems.push_back(el); + } + eat_blanks(); + } + } + + void eat_blanks() { + while (m_line.size()) { + if (m_line[0] == ' ') + m_line = m_line.substr(1); + else + break; + } + } + + void fill_formula_elem() { + fill_lisp_elem(m_formula_lisp_elem); + } + + void parse_line() { + if (m_line.find(":formula") == 0) { + int first_br = static_cast(m_line.find('(')); + if (first_br == -1) { + std::cout << "empty formula" << std::endl; + return; + } + m_line = m_line.substr(first_br); + fill_formula_elem(); + } + } + + void set_constraint_kind(formula_constraint & c, lisp_elem & el) { + if (el.m_head == "=") { + c.m_kind = EQ; + } else if (el.m_head == ">=") { + c.m_kind = GE; + } else if (el.m_head == "<=") { + c.m_kind = LE; + } else if (el.m_head == ">") { + c.m_kind = GT; + } else if (el.m_head == "<") { + c.m_kind = LT; + } else { + std::cout << "kind " << el.m_head << " is not supported " << std::endl; + set_error(); + } + } + + void adjust_rigth_side(formula_constraint & /* c*/, lisp_elem & /*el*/) { + // lean_assert(el.m_head == "0"); // do nothing for the time being + } + + void set_constraint_coeffs(formula_constraint & c, lisp_elem & el) { + lean_assert(el.m_elems.size() == 2); + set_constraint_coeffs_on_coeff_element(c, el.m_elems[0]); + adjust_rigth_side(c, el.m_elems[1]); + } + + + bool is_integer(std::string & s) { + if (s.size() == 0) return false; + return atoi(s.c_str()) != 0 || isdigit(s.c_str()[0]); + } + + void add_complex_sum_elem(formula_constraint & c, lisp_elem & el) { + if (el.m_head == "*") { + add_mult_elem(c, el.m_elems); + } else if (el.m_head == "~") { + lisp_elem & minel = el.m_elems[0]; + lean_assert(minel.is_simple()); + c.m_right_side += mpq(str_to_int(minel.m_head)); + } else { + std::cout << "unexpected input " << el.m_head << std::endl; + set_error(); + return; + } + } + + std::string get_name(lisp_elem & name) { + lean_assert(name.is_simple()); + lean_assert(!is_integer(name.m_head)); + return name.m_head; + } + + + void add_mult_elem(formula_constraint & c, std::vector & els) { + lean_assert(els.size() == 2); + mpq coeff = get_coeff(els[0]); + std::string col_name = get_name(els[1]); + c.add_pair(coeff, col_name); + } + + mpq get_coeff(lisp_elem & le) { + if (le.is_simple()) { + return mpq(str_to_int(le.m_head)); + } else { + lean_assert(le.m_head == "~"); + lean_assert(le.size() == 1); + lisp_elem & el = le.m_elems[0]; + lean_assert(el.is_simple()); + return -mpq(str_to_int(el.m_head)); + } + } + + int str_to_int(std::string & s) { + lean_assert(is_integer(s)); + return atoi(s.c_str()); + } + + void add_sum_elem(formula_constraint & c, lisp_elem & el) { + if (el.size()) { + add_complex_sum_elem(c, el); + } else { + lean_assert(is_integer(el.m_head)); + int v = atoi(el.m_head.c_str()); + mpq vr(v); + c.m_right_side -= vr; + } + } + + void add_sum(formula_constraint & c, std::vector & sum_els) { + for (auto & el : sum_els) + add_sum_elem(c, el); + } + + void set_constraint_coeffs_on_coeff_element(formula_constraint & c, lisp_elem & el) { + if (el.m_head == "*") { + add_mult_elem(c, el.m_elems); + } else if (el.m_head == "+") { + add_sum(c, el.m_elems); + } else { + lean_assert(false); // unexpected input + } + } + + void create_constraint(lisp_elem & el) { + formula_constraint c; + set_constraint_kind(c, el); + set_constraint_coeffs(c, el); + m_constraints.push_back(c); + } + + void fill_constraints() { + if (m_formula_lisp_elem.m_head != "and") { + std::cout << "unexpected top element " << m_formula_lisp_elem.m_head << std::endl; + set_error(); + return; + } + for (auto & el : m_formula_lisp_elem.m_elems) + create_constraint(el); + } + + void read() { + if (!m_file_stream.is_open()){ + std::cout << "cannot open file " << m_file_name << std::endl; + set_error(); + return; + } + while (m_is_OK && getline(m_file_stream, m_line)) { + parse_line(); + m_line_number++; + } + + m_file_stream.close(); + fill_constraints(); + } + + /* + void fill_lar_solver_on_row(row * row, lar_solver *solver) { + if (row->m_name != m_cost_row_name) { + lar_constraint c(get_lar_relation_from_row(row->m_type), row->m_right_side); + for (auto s : row->m_row_columns) { + var_index i = solver->add_var(s.first); + c.add_variable_to_constraint(i, s.second); + } + solver->add_constraint(&c); + } else { + // ignore the cost row + } + } + + + void fill_lar_solver_on_rows(lar_solver * solver) { + for (auto row_it : m_rows) { + fill_lar_solver_on_row(row_it.second, solver); + } + } + + void create_low_constraint_for_var(column* col, bound * b, lar_solver *solver) { + lar_constraint c(GE, b->m_low); + var_index i = solver->add_var(col->m_name); + c.add_variable_to_constraint(i, numeric_traits::one()); + solver->add_constraint(&c); + } + + void create_upper_constraint_for_var(column* col, bound * b, lar_solver *solver) { + lar_constraint c(LE, b->m_upper); + var_index i = solver->add_var(col->m_name); + c.add_variable_to_constraint(i, numeric_traits::one()); + solver->add_constraint(&c); + } + + void create_equality_contraint_for_var(column* col, bound * b, lar_solver *solver) { + lar_constraint c(EQ, b->m_fixed_value); + var_index i = solver->add_var(col->m_name); + c.add_variable_to_constraint(i, numeric_traits::one()); + solver->add_constraint(&c); + } + + void fill_lar_solver_on_columns(lar_solver * solver) { + for (auto s : m_columns) { + mps_reader::column * col = s.second; + solver->add_var(col->m_name); + auto b = col->m_bound; + if (b == nullptr) return; + + if (b->m_free) continue; + + if (b->m_low_is_set) { + create_low_constraint_for_var(col, b, solver); + } + if (b->m_upper_is_set) { + create_upper_constraint_for_var(col, b, solver); + } + if (b->m_value_is_fixed) { + create_equality_contraint_for_var(col, b, solver); + } + } + } + */ + + unsigned register_name(std::string s) { + auto it = m_name_to_var_index.find(s); + if (it!= m_name_to_var_index.end()) + return it->second; + + unsigned ret= m_name_to_var_index.size(); + m_name_to_var_index[s] = ret; + return ret; + } + + void add_constraint_to_solver(lar_solver * solver, formula_constraint & fc) { + vector> ls; + for (auto & it : fc.m_coeffs) { + ls.push_back(std::make_pair(it.first, solver->add_var(register_name(it.second)))); + } + solver->add_constraint(ls, fc.m_kind, fc.m_right_side); + } + + void fill_lar_solver(lar_solver * solver) { + for (formula_constraint & fc : m_constraints) + add_constraint_to_solver(solver, fc); + } + + + lar_solver * create_lar_solver() { + lar_solver * ls = new lar_solver(); + fill_lar_solver(ls); + return ls; + } + }; +} diff --git a/src/test/test_file_reader.h b/src/test/test_file_reader.h new file mode 100644 index 000000000..c7a9e3b8b --- /dev/null +++ b/src/test/test_file_reader.h @@ -0,0 +1,73 @@ +/* +Copyright (c) 2013 Microsoft Corporation. All rights reserved. +Released under Apache 2.0 license as described in the file LICENSE. + +Author: Lev Nachmanson +*/ +#pragma once + +// reads a text file +#include +#include +#include +#include +#include +#include "util/lp/lp_utils.h" +#include "util/lp/lp_solver.h" + +namespace lean { + +template +struct test_result { + lp_status m_status; + T m_cost; + std::unordered_map column_values; +}; + +template +class test_file_reader { + struct raw_blob { + std::vector m_unparsed_strings; + std::vector m_blobs; + }; + + struct test_file_blob { + std::string m_name; + std::string m_content; + std::unordered_map m_table; + std::unordered_map m_blobs; + + test_result * get_test_result() { + test_result * tr = new test_result(); + throw "not impl"; + return tr; + } + }; + std::ifstream m_file_stream; +public: + // constructor + test_file_reader(std::string file_name) : m_file_stream(file_name) { + if (!m_file_stream.is_open()) { + std::cout << "cannot open file " << "\'" << file_name << "\'" << std::endl; + } + } + + raw_blob scan_to_row_blob() { + } + + test_file_blob scan_row_blob_to_test_file_blob(raw_blob /* rblob */) { + } + + test_result * get_test_result() { + if (!m_file_stream.is_open()) { + return nullptr; + } + + raw_blob rblob = scan_to_row_blob(); + + test_file_blob tblob = scan_row_blob_to_test_file_blob(rblob); + + return tblob.get_test_result(); + } +}; +} diff --git a/src/util/lp/lp_params.pyg b/src/util/lp/lp_params.pyg new file mode 100644 index 000000000..3d31d0bdf --- /dev/null +++ b/src/util/lp/lp_params.pyg @@ -0,0 +1,13 @@ +def_module_params('lp', + export=True, + 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') + )) + + +