diff --git a/src/math/simplex/simplex_def.h b/src/math/simplex/simplex_def.h index 93b061dda..a3a44e34b 100644 --- a/src/math/simplex/simplex_def.h +++ b/src/math/simplex/simplex_def.h @@ -608,6 +608,7 @@ namespace simplex { // optimal return l_true; } + TRACE("simplex", tout << "x_i: v" << x_i << " x_j: v" << x_j << "\n";); var_info& vj = m_vars[x_j]; if (x_i == null_var) { if (inc && vj.m_upper_valid) { @@ -635,8 +636,9 @@ namespace simplex { // pivot(x_i, x_j, a_ij); - move_to_bound(x_i, inc == m.is_pos(a_ij)); + move_to_bound(x_i, inc != m.is_pos(a_ij)); SASSERT(well_formed_row(row(m_vars[x_j].m_base2row))); + TRACE("simplex", display(tout);); } return l_true; } @@ -651,22 +653,37 @@ namespace simplex { else { em.sub(vi.m_upper, vi.m_value, delta); } + TRACE("simplex", tout << "move " << (to_lower?"to_lower":"to_upper") + << " v" << x << " " << em.to_string(delta) << "\n";); col_iterator it = M.col_begin(x), end = M.col_end(x); for (; it != end && is_pos(delta); ++it) { + // + // base_coeff*s + coeff*x + R = 0 + // + // to_lower coeff > 0 base_coeff > 0 bound(s) + // ------------------------------------------------------ + // T T T !to_lower + // T T F to_lower + // T F T to_lower + // T F F !to_lower + // var_t s = m_row2base[it.get_row().id()]; var_info& vs = m_vars[s]; numeral const& coeff = it.get_row_entry().m_coeff; + numeral const& base_coeff = vs.m_base_coeff; SASSERT(!m.is_zero(coeff)); - bool inc_s = (m.is_pos(coeff) == to_lower); + bool base_to_lower = (m.is_pos(coeff) != m.is_pos(base_coeff)) == to_lower; eps_numeral const* bound = 0; - if (inc_s && vs.m_upper_valid) { + if (!base_to_lower && vs.m_upper_valid) { bound = &vs.m_upper; } - else if (!inc_s && vs.m_lower_valid) { + else if (base_to_lower && vs.m_lower_valid) { bound = &vs.m_lower; } if (bound) { + // |delta2*coeff| = |(bound-value)*base_coeff| em.sub(*bound, vs.m_value, delta2); + em.mul(delta2, base_coeff, delta2); em.div(delta2, coeff, delta2); abs(delta2); if (delta2 < delta) { @@ -694,8 +711,11 @@ namespace simplex { for (; it != end; ++it) { var_t x = it->m_var; if (x == v) continue; - bool is_pos = m.is_pos(it->m_coeff); + bool is_pos = m.is_pos(it->m_coeff) == m.is_pos(m_vars[v].m_base_coeff); if ((is_pos && at_upper(x)) || (!is_pos && at_lower(x))) { + TRACE("simplex", tout << "v" << x << " pos: " << is_pos + << " at upper: " << at_upper(x) + << " at lower: " << at_lower(x) << "\n";); continue; // variable cannot be used for improving bounds. // TBD check? } @@ -737,7 +757,11 @@ namespace simplex { var_info& vi = m_vars[s]; numeral const& a_ij = it.get_row_entry().m_coeff; numeral const& a_ii = vi.m_base_coeff; - bool inc_s = m.is_neg(a_ij) ? inc : !inc; + bool inc_s = (m.is_pos(a_ii) != m.is_pos(a_ij)) ? inc : !inc; + TRACE("simplex", tout << "v" << x_j << " incs: " << inc_s + << " upper valid:" << vi.m_upper_valid + << " lower valid:" << vi.m_lower_valid << "\n"; + display_row(tout, r);); if ((inc_s && !vi.m_upper_valid) || (!inc_s && !vi.m_lower_valid)) { continue; } @@ -751,14 +775,15 @@ namespace simplex { if (is_neg(curr_gain)) { curr_gain.neg(); } - if (x_i == null_var || (curr_gain < gain) || + if (x_i == null_var || (gain < curr_gain) || (is_zero(gain) && is_zero(curr_gain) && s < x_i)) { x_i = s; gain = curr_gain; new_a_ij = a_ij; + TRACE("simplex", tout << "x_j v" << x_j << " x_i v" << x_i << " gain: "; + tout << em.to_string(curr_gain) << "\n";); } } - TRACE("simplex", tout << "x_i v" << x_i << "\n";); return x_i; } diff --git a/src/smt/theory_diff_logic.h b/src/smt/theory_diff_logic.h index dff1b718c..70a3268d7 100644 --- a/src/smt/theory_diff_logic.h +++ b/src/smt/theory_diff_logic.h @@ -158,12 +158,13 @@ namespace smt { unsigned m_asserted_atoms_lim; unsigned m_asserted_qhead_old; }; + typedef dl_graph Graph; smt_params & m_params; arith_util m_util; arith_eq_adapter m_arith_eq_adapter; theory_diff_logic_statistics m_stats; - dl_graph m_graph; + Graph m_graph; theory_var m_zero; // cache the variable representing the zero variable. int_vector m_scc_id; // Cheap equality propagation eq_prop_info_set m_eq_prop_info_set; // set of existing equality prop infos @@ -356,9 +357,7 @@ namespace smt { void get_implied_bound_antecedents(edge_id bridge_edge, edge_id subsumed_edge, conflict_resolution & cr); - theory_var get_zero(sort* s) const { return m_zero; } - - theory_var get_zero(expr* e) const { return m_zero; } + theory_var get_zero() const { return m_zero; } void inc_conflicts(); diff --git a/src/smt/theory_diff_logic_def.h b/src/smt/theory_diff_logic_def.h index e0577f45a..553a4aef2 100644 --- a/src/smt/theory_diff_logic_def.h +++ b/src/smt/theory_diff_logic_def.h @@ -208,7 +208,7 @@ bool theory_diff_logic::internalize_atom(app * n, bool gate_ctx) { } else { target = mk_var(lhs); - source = get_zero(lhs); + source = get_zero(); } if (is_ge) { std::swap(target, source); @@ -698,7 +698,7 @@ theory_var theory_diff_logic::mk_num(app* n, rational const& r) { enode* e = 0; context& ctx = get_context(); if (r.is_zero()) { - v = get_zero(n); + v = get_zero(); } else if (ctx.e_internalized(n)) { e = ctx.get_enode(n); @@ -706,7 +706,7 @@ theory_var theory_diff_logic::mk_num(app* n, rational const& r) { SASSERT(v != null_theory_var); } else { - theory_var zero = get_zero(n); + theory_var zero = get_zero(); e = ctx.mk_enode(n, false, false, true); v = mk_var(e); // internalizer is marking enodes as interpreted whenever the associated ast is a value and a constant. @@ -1006,71 +1006,82 @@ inf_eps_rational theory_diff_logic::maximize(theory_var v) { IF_VERBOSE(1, for (unsigned i = 0; i < objective.size(); ++i) { - verbose_stream() << "Coefficient " << objective[i].second << " of theory_var " << objective[i].first << "\n"; + verbose_stream() << "Coefficient " << objective[i].second + << " of theory_var " << objective[i].first << "\n"; } verbose_stream() << "Free coefficient " << m_objective_consts[v] << "\n";); - -#if 0 - // disabled until fixed. - - - // Objective coefficients now become balances - vector balances(m_graph.get_num_nodes()); - balances.fill(fin_numeral::zero()); - fin_numeral sum = fin_numeral::zero(); - for (unsigned i = 0; i < objective.size(); ++i) { - fin_numeral balance(objective[i].second); - balances[objective[i].first] = balance; - sum += balance; - } - // HACK: assume that v0 is always value 0 - balances[0] = -sum; - - TRACE("arith", display(tout);); - network_flow net_flow(m_graph, balances); - min_flow_result result = net_flow.min_cost(); - SASSERT(result != UNBOUNDED); - if (result == OPTIMAL) { - numeral objective_value = net_flow.get_optimal_solution(m_objective_assignments[v], true) + numeral(m_objective_consts[v]); - IF_VERBOSE(1, verbose_stream() << "Optimal value of objective " << v << ": " << objective_value << std::endl;); - - DEBUG_CODE( - numeral initial_value = numeral(m_objective_consts[v]); - for (unsigned i = 0; i < objective.size(); ++i) { - initial_value += fin_numeral(objective[i].second) * m_graph.get_assignment(objective[i].first); - } - IF_VERBOSE(1, verbose_stream() << "Initial value of objective " << v << ": " << initial_value << std::endl;); - // FIXME: Network Simplex lose precisions when handling infinitesimals - SASSERT(objective_value >= initial_value.get_rational());); - vector & current_assigments = m_objective_assignments[v]; - SASSERT(!current_assigments.empty()); - // Normalize optimal assignments so that v0 is fixed to 0 - for (unsigned i = 1; i < current_assigments.size(); ++i) { - current_assigments[i] -= current_assigments[0]; - } - - ast_manager& m = get_manager(); - IF_VERBOSE(1, - verbose_stream() << "Optimal assigment:" << std::endl; - for (unsigned i = 0; i < objective.size(); ++i) { - theory_var v = objective[i].first; - verbose_stream() << mk_pp(get_enode(v)->get_owner(), m) << " |-> " << current_assigments[v] << std::endl; - }); - rational r = objective_value.get_rational().to_rational(); - rational i = objective_value.get_infinitesimal().to_rational(); - return inf_eps_rational(inf_rational(r, i)); - } - else { - // Dual problem is infeasible, primal problem is unbounded - SASSERT(result == INFEASIBLE); - IF_VERBOSE(1, verbose_stream() << "Unbounded objective" << std::endl;); - return inf_eps_rational::infinity(); + unsigned num_nodes = m_graph.get_num_nodes(); + unsigned num_edges = m_graph.get_num_edges(); + vector > const& es = m_graph.get_all_edges(); + S.ensure_var(num_nodes + num_edges + m_objectives.size()); + for (unsigned i = 0; i < num_nodes; ++i) { + numeral const& a = m_graph.get_assignment(i); + rational fin = a.get_rational().to_rational(); + rational inf = a.get_infinitesimal().to_rational(); + mpq_inf q(fin.to_mpq(), inf.to_mpq()); + S.set_value(i, q); } + S.set_lower(get_zero(), mpq_inf(mpq(0), mpq(0))); + S.set_upper(get_zero(), mpq_inf(mpq(0), mpq(0))); + svector vars; + unsynch_mpq_manager mgr; + scoped_mpq_vector coeffs(mgr); + coeffs.push_back(mpq(1)); + coeffs.push_back(mpq(-1)); + coeffs.push_back(mpq(-1)); + vars.resize(3); + for (unsigned i = 0; i < es.size(); ++i) { + dl_edge const& e = es[i]; + if (e.is_enabled()) { + unsigned base_var = num_nodes + i; + vars[0] = e.get_target(); + vars[1] = e.get_source(); + vars[2] = base_var; + S.add_row(base_var, 3, vars.c_ptr(), coeffs.c_ptr()); + // t - s <= w + // t - s - b = 0, b >= w + numeral const& w = e.get_weight(); + rational fin = w.get_rational().to_rational(); + rational inf = w.get_infinitesimal().to_rational(); + mpq_inf q(fin.to_mpq(),inf.to_mpq()); + S.set_upper(base_var, q); + } + } + unsigned w = num_nodes + num_edges + v; -#endif - return inf_eps_rational::infinity(); + // add objective function as row. + coeffs.reset(); + vars.reset(); + for (unsigned i = 0; i < objective.size(); ++i) { + coeffs.push_back(objective[i].second.to_mpq()); + vars.push_back(objective[i].first); + } + coeffs.push_back(mpq(1)); + vars.push_back(w); + S.add_row(w, vars.size(), vars.c_ptr(), coeffs.c_ptr()); + TRACE("opt", S.display(tout); display(tout);); + + // optimize + lbool is_sat = S.make_feasible(); + if (is_sat == l_undef) { + return inf_eps_rational::infinity(); + } + TRACE("opt", S.display(tout); ); + SASSERT(is_sat != l_false); + lbool is_fin = S.minimize(w); + switch (is_fin) { + case l_true: { + simplex::mpq_ext::eps_numeral const& val = S.get_value(w); + inf_rational r(-rational(val.first), -rational(val.second)); + TRACE("opt", tout << r << " " << "\n"; ); + return inf_eps_rational(rational(0), r); + } + default: + TRACE("opt", tout << "unbounded\n"; ); + return inf_eps_rational::infinity(); + } } template