diff --git a/src/api/api_opt.cpp b/src/api/api_opt.cpp
index afc33e396..b1640506e 100644
--- a/src/api/api_opt.cpp
+++ b/src/api/api_opt.cpp
@@ -20,6 +20,7 @@ Revision History:
#include"api_log_macros.h"
#include"api_context.h"
#include"api_util.h"
+#include"api_model.h"
#include"opt_context.h"
#include"cancel_eh.h"
#include"scoped_timer.h"
@@ -122,6 +123,23 @@ extern "C" {
Z3_CATCH_RETURN(Z3_L_UNDEF);
}
+ Z3_model Z3_API Z3_optimize_get_model(Z3_context c, Z3_optimize o) {
+ Z3_TRY;
+ LOG_Z3_optimize_get_model(c, o);
+ RESET_ERROR_CODE();
+ model_ref _m;
+ to_optimize_ref(o).get_model(_m);
+ if (!_m) {
+ SET_ERROR_CODE(Z3_INVALID_USAGE);
+ RETURN_Z3(0);
+ }
+ Z3_model_ref * m_ref = alloc(Z3_model_ref);
+ m_ref->m_model = _m;
+ mk_c(c)->save_object(m_ref);
+ RETURN_Z3(of_model(m_ref));
+ Z3_CATCH_RETURN(0);
+ }
+
void Z3_API Z3_optimize_set_params(Z3_context c, Z3_optimize o, Z3_params p) {
Z3_TRY;
LOG_Z3_optimize_set_params(c, o, p);
diff --git a/src/api/dotnet/Optimize.cs b/src/api/dotnet/Optimize.cs
index 3b54dca20..4526ead88 100644
--- a/src/api/dotnet/Optimize.cs
+++ b/src/api/dotnet/Optimize.cs
@@ -108,6 +108,26 @@ namespace Microsoft.Z3
default: return Status.UNKNOWN;
}
}
+
+ ///
+ /// The model of the last Check.
+ ///
+ ///
+ /// The result is null if Check was not invoked before,
+ /// if its results was not SATISFIABLE, or if model production is not enabled.
+ ///
+ public Model Model
+ {
+ get
+ {
+ IntPtr x = Native.Z3_optimize_get_model(Context.nCtx, NativeObject);
+ if (x == IntPtr.Zero)
+ return null;
+ else
+ return new Model(Context, x);
+ }
+ }
+
public void Maximize(ArithExpr e) {
Native.Z3_optimize_maximize(Context.nCtx, NativeObject, e.NativeObject);
diff --git a/src/api/z3_api.h b/src/api/z3_api.h
index f6064d543..a52ad6106 100644
--- a/src/api/z3_api.h
+++ b/src/api/z3_api.h
@@ -6021,6 +6021,17 @@ END_MLAPI_EXCLUDE
Z3_lbool Z3_API Z3_optimize_check(Z3_context c, Z3_optimize o);
+ /**
+ \brief Retrieve the model for the last #Z3_optimize_check
+
+ The error handler is invoked if a model is not available because
+ the commands above were not invoked for the given optimization
+ solver, or if the result was \c Z3_L_FALSE.
+
+ def_API('Z3_optimize_get_model', MODEL, (_in(CONTEXT), _in(OPTIMIZE)))
+ */
+ Z3_model Z3_API Z3_optimize_get_model(Z3_context c, Z3_optimize o);
+
/**
\brief Set parameters on optimization context.
diff --git a/src/opt/opt_context.cpp b/src/opt/opt_context.cpp
index fafd88ecf..31202927c 100644
--- a/src/opt/opt_context.cpp
+++ b/src/opt/opt_context.cpp
@@ -82,6 +82,10 @@ namespace opt {
}
}
+ void context::get_model(model_ref& mdl) {
+ get_solver().get_model(mdl);
+ }
+
lbool context::execute_min_max(unsigned index, bool committed, bool is_max) {
// HACK: reuse m_optsmt without regard for box reuse and not considering
// use-case of lex.
diff --git a/src/opt/opt_context.h b/src/opt/opt_context.h
index 9a0b49541..a1f391d1d 100644
--- a/src/opt/opt_context.h
+++ b/src/opt/opt_context.h
@@ -75,6 +75,8 @@ namespace opt {
lbool optimize();
+ void get_model(model_ref& m);
+
void set_cancel(bool f);
void reset_cancel() { set_cancel(false); }
void cancel() { set_cancel(true); }
diff --git a/src/smt/network_flow.h b/src/smt/network_flow.h
index 0526de965..65575965f 100644
--- a/src/smt/network_flow.h
+++ b/src/smt/network_flow.h
@@ -77,9 +77,9 @@ namespace smt {
protected:
graph & m_graph;
svector & m_states;
- vector & m_potentials;
- edge_id & m_enter_id;
-
+ vector & m_potentials;
+ edge_id & m_enter_id;
+ bool edge_in_tree(edge_id id) const { return m_states[id] == BASIS; }
public:
pivot_rule_impl(graph & g, vector & potentials,
svector & states, edge_id & enter_id)
@@ -88,46 +88,21 @@ namespace smt {
m_states(states),
m_enter_id(enter_id) {
}
+ virtual ~pivot_rule_impl() {}
virtual bool choose_entering_edge() = 0;
+ virtual pivot_rule rule() const = 0;
};
class first_eligible_pivot : public pivot_rule_impl {
- private:
edge_id m_next_edge;
-
public:
first_eligible_pivot(graph & g, vector & potentials,
svector & states, edge_id & enter_id) :
pivot_rule_impl(g, potentials, states, enter_id),
m_next_edge(0) {
}
-
- bool choose_entering_edge() {
- TRACE("network_flow", tout << "choose_entering_edge...\n";);
- int num_edges = m_graph.get_num_edges();
- for (int i = m_next_edge; i < m_next_edge + num_edges; ++i) {
- edge_id id = (i >= num_edges) ? (i - num_edges) : i;
- node src = m_graph.get_source(id);
- node tgt = m_graph.get_target(id);
- if (m_states[id] != BASIS) {
- numeral cost = m_potentials[src] - m_potentials[tgt] - m_graph.get_weight(id);
- if (cost.is_pos()) {
- m_enter_id = id;
- TRACE("network_flow", {
- tout << "Found entering edge " << id << " between node ";
- tout << src << " and node " << tgt << " with reduced cost = " << cost << "...\n";
- });
- m_next_edge = m_enter_id;
- if (m_next_edge >= num_edges) {
- m_next_edge -= num_edges;
- }
- return true;
- }
- }
- }
- TRACE("network_flow", tout << "Found no entering edge...\n";);
- return false;
- };
+ virtual bool choose_entering_edge();
+ virtual pivot_rule rule() const { return FIRST_ELIGIBLE; }
};
class best_eligible_pivot : public pivot_rule_impl {
@@ -136,33 +111,8 @@ namespace smt {
svector & states, edge_id & enter_id) :
pivot_rule_impl(g, potentials, states, enter_id) {
}
-
- bool choose_entering_edge() {
- TRACE("network_flow", tout << "choose_entering_edge...\n";);
- unsigned num_edges = m_graph.get_num_edges();
- numeral max = numeral::zero();
- for (unsigned i = 0; i < num_edges; ++i) {
- node src = m_graph.get_source(i);
- node tgt = m_graph.get_target(i);
- if (m_states[i] != BASIS) {
- numeral cost = m_potentials[src] - m_potentials[tgt] - m_graph.get_weight(i);
- if (cost > max) {
- max = cost;
- m_enter_id = i;
- }
- }
- }
- if (max.is_pos()) {
- TRACE("network_flow", {
- tout << "Found entering edge " << m_enter_id << " between node ";
- tout << m_graph.get_source(m_enter_id) << " and node " << m_graph.get_target(m_enter_id);
- tout << " with reduced cost = " << max << "...\n";
- });
- return true;
- }
- TRACE("network_flow", tout << "Found no entering edge...\n";);
- return false;
- };
+ virtual pivot_rule rule() const { return BEST_ELIGIBLE; }
+ virtual bool choose_entering_edge();
};
class candidate_list_pivot : public pivot_rule_impl {
@@ -186,94 +136,22 @@ namespace smt {
m_candidates(m_num_candidates) {
}
- bool choose_entering_edge() {
- if (m_current_length == 0 || m_minor_step == MINOR_STEP_LIMIT) {
- // Build the candidate list
- unsigned num_edges = m_graph.get_num_edges();
- numeral max = numeral::zero();
- m_current_length = 0;
- for (unsigned i = m_next_edge; i < m_next_edge + num_edges; ++i) {
- edge_id id = (i >= num_edges) ? i - num_edges : i;
- node src = m_graph.get_source(id);
- node tgt = m_graph.get_target(id);
- if (m_states[id] != BASIS) {
- numeral cost = m_potentials[src] - m_potentials[tgt] - m_graph.get_weight(id);
- if (cost.is_pos()) {
- m_candidates[m_current_length] = id;
- ++m_current_length;
- if (cost > max) {
- max = cost;
- m_enter_id = id;
- }
- }
- if (m_current_length >= m_num_candidates) break;
- }
- }
- m_next_edge = m_enter_id;
- m_minor_step = 1;
- if (max.is_pos()) {
- TRACE("network_flow", {
- tout << "Found entering edge " << m_enter_id << " between node ";
- tout << m_graph.get_source(m_enter_id) << " and node " << m_graph.get_target(m_enter_id);
- tout << " with reduced cost = " << max << "...\n";
- });
- return true;
- }
- TRACE("network_flow", tout << "Found no entering edge...\n";);
- return false;
- }
+ virtual pivot_rule rule() const { return CANDIDATE_LIST; }
- ++m_minor_step;
- numeral max = numeral::zero();
- for (unsigned i = 0; i < m_current_length; ++i) {
- edge_id id = m_candidates[i];
- node src = m_graph.get_source(id);
- node tgt = m_graph.get_target(id);
- if (m_states[id] != BASIS) {
- numeral cost = m_potentials[src] - m_potentials[tgt] - m_graph.get_weight(id);
- if (cost > max) {
- max = cost;
- m_enter_id = id;
- }
- // Remove stale candidates
- if (!cost.is_pos()) {
- --m_current_length;
- m_candidates[i] = m_candidates[m_current_length];
- --i;
- }
- }
- }
- if (max.is_pos()) {
- TRACE("network_flow", {
- tout << "Found entering edge " << m_enter_id << " between node ";
- tout << m_graph.get_source(m_enter_id) << " and node " << m_graph.get_target(m_enter_id);
- tout << " with reduced cost = " << max << "...\n";
- });
- return true;
- }
- TRACE("network_flow", tout << "Found no entering edge...\n";);
- return false;
- };
+ virtual bool choose_entering_edge();
};
- graph m_graph;
- spanning_tree_base * m_tree;
-
- // Denote supply/demand b_i on node i
- vector m_balances;
-
- // Duals of flows which are convenient to compute dual solutions
- vector m_potentials;
-
- // Basic feasible flows
- vector m_flows;
-
- svector m_states;
-
- unsigned m_step;
-
- edge_id m_enter_id, m_leave_id;
- optional m_delta;
+ graph m_graph;
+ scoped_ptr m_tree;
+ scoped_ptr m_pivot;
+ vector m_balances; // Denote supply/demand b_i on node i
+ vector m_potentials; // Duals of flows which are convenient to compute dual solutions
+ vector m_flows; // Basic feasible flows
+ svector m_states;
+ unsigned m_step;
+ edge_id m_enter_id;
+ edge_id m_leave_id;
+ optional m_delta;
// Initialize the network with a feasible spanning tree
void initialize();
@@ -290,6 +168,8 @@ namespace smt {
void update_spanning_tree();
+ numeral get_cost() const;
+
bool edge_in_tree(edge_id id) const;
bool is_infeasible();
diff --git a/src/smt/network_flow_def.h b/src/smt/network_flow_def.h
index 72fedbac9..683ea0b98 100644
--- a/src/smt/network_flow_def.h
+++ b/src/smt/network_flow_def.h
@@ -26,6 +26,128 @@ Notes:
namespace smt {
+ template
+ bool network_flow::first_eligible_pivot::choose_entering_edge() {
+ TRACE("network_flow", tout << "choose_entering_edge...\n";);
+ int num_edges = m_graph.get_num_edges();
+ for (int i = m_next_edge; i < m_next_edge + num_edges; ++i) {
+ edge_id id = i % num_edges;
+ if (edge_in_tree(id)) {
+ continue;
+ }
+ node src = m_graph.get_source(id);
+ node tgt = m_graph.get_target(id);
+ numeral cost = m_potentials[src] - m_potentials[tgt] - m_graph.get_weight(id);
+ if (cost.is_pos()) {
+ m_next_edge = m_enter_id = id;
+ TRACE("network_flow",
+ tout << "Found entering edge " << id << " between node ";
+ tout << src << " and node " << tgt << " with reduced cost = " << cost << "...\n";);
+ return true;
+ }
+ }
+ TRACE("network_flow", tout << "Found no entering edge...\n";);
+ return false;
+ };
+
+ template
+ bool network_flow::best_eligible_pivot::choose_entering_edge() {
+ TRACE("network_flow", tout << "choose_entering_edge...\n";);
+ unsigned num_edges = m_graph.get_num_edges();
+ numeral max = numeral::zero();
+ for (unsigned i = 0; i < num_edges; ++i) {
+ node src = m_graph.get_source(i);
+ node tgt = m_graph.get_target(i);
+ if (!edge_in_tree(i)) {
+ numeral cost = m_potentials[src] - m_potentials[tgt] - m_graph.get_weight(i);
+ if (cost > max) {
+ max = cost;
+ m_enter_id = i;
+ }
+ }
+ }
+ if (max.is_pos()) {
+ TRACE("network_flow", {
+ tout << "Found entering edge " << m_enter_id << " between node ";
+ tout << m_graph.get_source(m_enter_id) << " and node " << m_graph.get_target(m_enter_id);
+ tout << " with reduced cost = " << max << "...\n";
+ });
+ return true;
+ }
+ TRACE("network_flow", tout << "Found no entering edge...\n";);
+ return false;
+ };
+
+ template
+ bool network_flow::candidate_list_pivot::choose_entering_edge() {
+ if (m_current_length == 0 || m_minor_step == MINOR_STEP_LIMIT) {
+ // Build the candidate list
+ unsigned num_edges = m_graph.get_num_edges();
+ numeral max = numeral::zero();
+ m_current_length = 0;
+ for (unsigned i = m_next_edge; i < m_next_edge + num_edges; ++i) {
+ edge_id id = (i >= num_edges) ? i - num_edges : i;
+ node src = m_graph.get_source(id);
+ node tgt = m_graph.get_target(id);
+ if (!edge_in_tree(id)) {
+ numeral cost = m_potentials[src] - m_potentials[tgt] - m_graph.get_weight(id);
+ if (cost.is_pos()) {
+ m_candidates[m_current_length] = id;
+ ++m_current_length;
+ if (cost > max) {
+ max = cost;
+ m_enter_id = id;
+ }
+ }
+ if (m_current_length >= m_num_candidates) break;
+ }
+ }
+ m_next_edge = m_enter_id;
+ m_minor_step = 1;
+ if (max.is_pos()) {
+ TRACE("network_flow", {
+ tout << "Found entering edge " << m_enter_id << " between node ";
+ tout << m_graph.get_source(m_enter_id) << " and node " << m_graph.get_target(m_enter_id);
+ tout << " with reduced cost = " << max << "...\n";
+ });
+ return true;
+ }
+ TRACE("network_flow", tout << "Found no entering edge...\n";);
+ return false;
+ }
+
+ ++m_minor_step;
+ numeral max = numeral::zero();
+ for (unsigned i = 0; i < m_current_length; ++i) {
+ edge_id id = m_candidates[i];
+ node src = m_graph.get_source(id);
+ node tgt = m_graph.get_target(id);
+ if (!edge_in_tree(id)) {
+ numeral cost = m_potentials[src] - m_potentials[tgt] - m_graph.get_weight(id);
+ if (cost > max) {
+ max = cost;
+ m_enter_id = id;
+ }
+ // Remove stale candidates
+ if (!cost.is_pos()) {
+ --m_current_length;
+ m_candidates[i] = m_candidates[m_current_length];
+ --i;
+ }
+ }
+ }
+ if (max.is_pos()) {
+ TRACE("network_flow", {
+ tout << "Found entering edge " << m_enter_id << " between node ";
+ tout << m_graph.get_source(m_enter_id) << " and node " << m_graph.get_target(m_enter_id);
+ tout << " with reduced cost = " << max << "...\n";
+ });
+ return true;
+ }
+ TRACE("network_flow", tout << "Found no entering edge...\n";);
+ return false;
+ };
+
template
network_flow::network_flow(graph & g, vector const & balances) :
m_balances(balances) {
@@ -42,7 +164,7 @@ namespace smt {
}
}
TRACE("network_flow", {
- tout << "Different logic optimization:" << std::endl;
+ tout << "Difference logic optimization:" << std::endl;
display_dual(tout);
tout << "Minimum cost flow:" << std::endl;
display(tout);
@@ -52,6 +174,276 @@ namespace smt {
m_tree = alloc(basic_spanning_tree, m_graph);
}
+
+ template
+ void network_flow::initialize() {
+ TRACE("network_flow", tout << "initialize...\n";);
+ // Create an artificial root node to construct initial spanning tree
+ unsigned num_nodes = m_graph.get_num_nodes();
+ unsigned num_edges = m_graph.get_num_edges();
+
+ node root = num_nodes;
+ m_graph.init_var(root);
+
+ m_potentials.resize(num_nodes + 1);
+ m_potentials[root] = numeral::zero();
+
+ m_balances.resize(num_nodes + 1);
+ fin_numeral sum_supply = fin_numeral::zero();
+ for (unsigned i = 0; i < num_nodes; ++i) {
+ sum_supply += m_balances[i];
+ }
+ m_balances[root] = -sum_supply;
+
+ m_flows.resize(num_nodes + num_edges);
+ m_flows.fill(numeral::zero());
+ m_states.resize(num_nodes + num_edges);
+ m_states.fill(LOWER);
+
+ // Create artificial edges from/to root node to/from other nodes and initialize the spanning tree
+ svector tree;
+ for (unsigned i = 0; i < num_nodes; ++i) {
+ bool is_forward = !m_balances[i].is_neg();
+ m_states[num_edges + i] = BASIS;
+ node src = is_forward ? i : root;
+ node tgt = is_forward ? root : i;
+ m_flows[num_edges + i] = is_forward ? m_balances[i] : -m_balances[i];
+ m_potentials[i] = is_forward ? numeral::one() : -numeral::one();
+ tree.push_back(m_graph.add_edge(src, tgt, numeral::one(), explanation()));
+ }
+
+ m_tree->initialize(tree);
+
+ TRACE("network_flow",
+ tout << pp_vector("Potentials", m_potentials);
+ tout << pp_vector("Flows", m_flows);
+ tout << "Cost: " << get_cost() << "\n";
+ tout << "Spanning tree:\n";
+ display_spanning_tree(tout););
+ SASSERT(check_well_formed());
+ }
+
+ template
+ void network_flow::update_potentials() {
+ node src = m_graph.get_source(m_enter_id);
+ node tgt = m_graph.get_target(m_enter_id);
+ numeral cost = m_potentials[src] - m_potentials[tgt] - m_graph.get_weight(m_enter_id);
+ numeral change;
+ node start;
+ if (m_tree->in_subtree_t2(tgt)) {
+ change = cost;
+ start = tgt;
+ }
+ else {
+ change = -cost;
+ start = src;
+ }
+ SASSERT(m_tree->in_subtree_t2(start));
+ TRACE("network_flow", tout << "update_potentials of T_" << start << " with change = " << change << "...\n";);
+ svector descendants;
+ m_tree->get_descendants(start, descendants);
+ SASSERT(descendants.size() >= 1);
+ for (unsigned i = 0; i < descendants.size(); ++i) {
+ node u = descendants[i];
+ m_potentials[u] += change;
+ }
+ TRACE("network_flow", tout << pp_vector("Potentials", m_potentials););
+ }
+
+ template
+ void network_flow::update_flows() {
+ TRACE("network_flow", tout << "update_flows...\n";);
+ m_flows[m_enter_id] += *m_delta;
+ node src = m_graph.get_source(m_enter_id);
+ node tgt = m_graph.get_target(m_enter_id);
+ svector path;
+ svector against;
+ m_tree->get_path(src, tgt, path, against);
+ SASSERT(path.size() >= 1);
+ for (unsigned i = 0; i < path.size(); ++i) {
+ edge_id e_id = path[i];
+ m_flows[e_id] += against[i] ? - *m_delta : *m_delta;
+ }
+ TRACE("network_flow", tout << pp_vector("Flows", m_flows););
+ }
+
+ template
+ bool network_flow::choose_leaving_edge() {
+ TRACE("network_flow", tout << "choose_leaving_edge...\n";);
+ node src = m_graph.get_source(m_enter_id);
+ node tgt = m_graph.get_target(m_enter_id);
+ m_delta.set_invalid();
+ edge_id leave_id = null_edge_id;
+ svector path;
+ svector against;
+ m_tree->get_path(src, tgt, path, against);
+ SASSERT(path.size() >= 1);
+ for (unsigned i = 0; i < path.size(); ++i) {
+ edge_id e_id = path[i];
+ if (against[i] && (!m_delta || m_flows[e_id] < *m_delta)) {
+ m_delta = m_flows[e_id];
+ leave_id = e_id;
+ }
+ }
+ m_leave_id = leave_id;
+
+ TRACE("network_flow",
+ tout << "Leaving edge " << m_leave_id;
+ tout << " between node " << m_graph.get_source(m_leave_id);
+ tout << " and node " << m_graph.get_target(m_leave_id) ;
+ if (m_delta) tout << " with delta = " << *m_delta;
+ tout << "\n";);
+
+ return m_delta;
+ }
+
+ template
+ void network_flow::update_spanning_tree() {
+ m_tree->update(m_enter_id, m_leave_id);
+ }
+
+ template
+ bool network_flow::choose_entering_edge(pivot_rule pr) {
+ if (!m_pivot || pr != m_pivot->rule()) {
+ switch (pr) {
+ case FIRST_ELIGIBLE:
+ m_pivot = alloc(first_eligible_pivot, m_graph, m_potentials, m_states, m_enter_id);
+ break;
+ case BEST_ELIGIBLE:
+ m_pivot = alloc(best_eligible_pivot, m_graph, m_potentials, m_states, m_enter_id);
+ break;
+ case CANDIDATE_LIST:
+ m_pivot = alloc(candidate_list_pivot, m_graph, m_potentials, m_states, m_enter_id);
+ break;
+ default:
+ UNREACHABLE();
+ }
+ }
+ return m_pivot->choose_entering_edge();
+ }
+
+ // Minimize cost flows
+ template
+ min_flow_result network_flow::min_cost(pivot_rule pr) {
+ initialize();
+ while (choose_entering_edge(pr)) {
+ bool bounded = choose_leaving_edge();
+ if (!bounded) return UNBOUNDED;
+ update_flows();
+ if (m_enter_id != m_leave_id) {
+ SASSERT(edge_in_tree(m_leave_id));
+ SASSERT(!edge_in_tree(m_enter_id));
+ m_states[m_enter_id] = BASIS;
+ m_states[m_leave_id] = LOWER;
+ update_spanning_tree();
+ update_potentials();
+ TRACE("network_flow",
+ tout << "Spanning tree:\n";
+ display_spanning_tree(tout);
+ tout << "Cost: " << get_cost() << "\n";);
+ SASSERT(check_well_formed());
+ }
+ }
+ if (is_infeasible()) return INFEASIBLE;
+ TRACE("network_flow", tout << "Found optimal solution.\n";);
+ SASSERT(check_optimal());
+ return OPTIMAL;
+ }
+
+ template
+ bool network_flow::is_infeasible() {
+#if 0
+ // Flows of artificial arcs should be zero
+ unsigned num_nodes = m_graph.get_num_nodes();
+ unsigned num_edges = m_graph.get_num_edges();
+ for (unsigned i = 1; i < num_nodes; ++i) {
+ if (m_flows[num_edges - i].is_pos()) return true;
+ }
+#endif
+ return false;
+ }
+
+ // Get the optimal solution
+ template
+ typename network_flow::numeral network_flow::get_optimal_solution(vector & result, bool is_dual) {
+ numeral objective_value = get_cost();
+ result.reset();
+ if (is_dual) {
+ result.append(m_potentials);
+ }
+ else {
+ result.append(m_flows);
+ }
+ return objective_value;
+ }
+
+ template
+ typename network_flow::numeral network_flow::get_cost() const {
+ numeral objective_value = numeral::zero();
+ unsigned num_edges = m_graph.get_num_edges();
+ for (unsigned i = 0; i < num_edges; ++i) {
+ if (edge_in_tree(i)) {
+ objective_value += m_flows[i].get_rational() * m_graph.get_weight(i);
+ }
+ }
+ return objective_value;
+ }
+
+ template
+ bool network_flow::edge_in_tree(edge_id id) const {
+ return m_states[id] == BASIS;
+ }
+
+ template
+ bool network_flow::check_well_formed() {
+ SASSERT(m_tree->check_well_formed());
+ SASSERT(!m_delta || !(*m_delta).is_neg());
+
+ // m_flows are zero on non-basic edges
+ for (unsigned i = 0; i < m_flows.size(); ++i) {
+ SASSERT(!m_flows[i].is_neg());
+ SASSERT(edge_in_tree(i) || m_flows[i].is_zero());
+ }
+
+ unsigned num_edges = m_graph.get_num_edges();
+ for (unsigned i = 0; i < num_edges; ++i) {
+ if (edge_in_tree(i)) {
+ dl_var src = m_graph.get_source(i);
+ dl_var tgt = m_graph.get_target(i);
+ numeral weight = m_graph.get_weight(i);
+ SASSERT(m_potentials[src] - m_potentials[tgt] == weight);
+ }
+ }
+
+ return true;
+ }
+
+ template
+ bool network_flow::check_optimal() {
+ numeral total_cost = get_cost();
+ unsigned num_edges = m_graph.get_num_edges();
+
+ for (unsigned i = 0; i < num_edges; ++i) {
+ dl_var src = m_graph.get_source(i);
+ dl_var tgt = m_graph.get_target(i);
+ numeral weight = m_graph.get_weight(i);
+ SASSERT(m_potentials[src] - m_potentials[tgt] <= weight);
+ }
+
+ // m_flows are zero on non-basic edges
+ for (unsigned i = 0; i < m_flows.size(); ++i) {
+ SASSERT(edge_in_tree(i) || m_flows[i].is_zero());
+ }
+ numeral total_balance = numeral::zero();
+ for (unsigned i = 0; i < m_potentials.size(); ++i) {
+ total_balance += m_balances[i] * m_potentials[i];
+ }
+ TRACE("network_flow", tout << "Total balance: " << total_balance << ", total cost: " << total_cost << std::endl;);
+ return total_cost == total_balance;
+ }
+
+ // display methods
+
template
void network_flow::display(std::ofstream & os) {
vector const & es = m_graph.get_all_edges();
@@ -96,6 +488,7 @@ namespace smt {
os << "(optimize)" << std::endl;
}
+
template
void network_flow::display_dual(std::ofstream & os) {
for (unsigned i = 0; i < m_graph.get_num_nodes(); ++i) {
@@ -115,296 +508,31 @@ namespace smt {
os << "(optimize)" << std::endl;
}
- template
- void network_flow::initialize() {
- TRACE("network_flow", tout << "initialize...\n";);
- // Create an artificial root node to construct initial spanning tree
- unsigned num_nodes = m_graph.get_num_nodes();
- unsigned num_edges = m_graph.get_num_edges();
-
- node root = num_nodes;
- m_graph.init_var(root);
-
- m_potentials.resize(num_nodes + 1);
- m_potentials[root] = numeral::zero();
-
- m_balances.resize(num_nodes + 1);
- fin_numeral sum_supply = fin_numeral::zero();
- for (unsigned i = 0; i < num_nodes; ++i) {
- sum_supply += m_balances[i];
- }
- m_balances[root] = -sum_supply;
-
- m_flows.resize(num_nodes + num_edges);
- m_flows.fill(numeral::zero());
- m_states.resize(num_nodes + num_edges);
- m_states.fill(LOWER);
-
- // Create artificial edges from/to root node to/from other nodes and initialize the spanning tree
- svector tree;
- bool is_forward;
- for (unsigned i = 0; i < num_nodes; ++i) {
- is_forward = !m_balances[i].is_neg();
- m_states[num_edges + i] = BASIS;
- node src = is_forward ? i : root;
- node tgt = is_forward ? root : i;
- m_flows[num_edges + i] = is_forward ? m_balances[i] : -m_balances[i];
- m_potentials[i] = is_forward ? numeral::one() : -numeral::one();
- tree.push_back(m_graph.add_edge(src, tgt, numeral::one(), explanation()));
- }
-
- m_tree->initialize(tree);
-
- TRACE("network_flow", {
- tout << pp_vector("Potentials", m_potentials, true) << pp_vector("Flows", m_flows);
- });
- TRACE("network_flow", tout << "Spanning tree:\n"; display_spanning_tree(tout););
- SASSERT(check_well_formed());
- }
-
- template
- void network_flow::update_potentials() {
- node src = m_graph.get_source(m_enter_id);
- node tgt = m_graph.get_target(m_enter_id);
- numeral cost = m_potentials[src] - m_potentials[tgt] - m_graph.get_weight(m_enter_id);
- numeral change;
- node start;
- if (m_tree->in_subtree_t2(tgt)) {
- change = cost;
- start = tgt;
- }
- else {
- change = -cost;
- start = src;
- }
- SASSERT(m_tree->in_subtree_t2(start));
- TRACE("network_flow", tout << "update_potentials of T_" << start << " with change = " << change << "...\n";);
- svector descendants;
- m_tree->get_descendants(start, descendants);
- SASSERT(descendants.size() >= 1);
- for (unsigned i = 0; i < descendants.size(); ++i) {
- node u = descendants[i];
- m_potentials[u] += change;
- }
- TRACE("network_flow", tout << pp_vector("Potentials", m_potentials, true););
- }
-
- template
- void network_flow::update_flows() {
- TRACE("network_flow", tout << "update_flows...\n";);
- m_flows[m_enter_id] += *m_delta;
- node src = m_graph.get_source(m_enter_id);
- node tgt = m_graph.get_target(m_enter_id);
- svector path;
- svector against;
- m_tree->get_path(src, tgt, path, against);
- SASSERT(path.size() >= 1);
- for (unsigned i = 0; i < path.size(); ++i) {
- edge_id e_id = path[i];
- m_flows[e_id] += against[i] ? - *m_delta : *m_delta;
- }
- TRACE("network_flow", tout << pp_vector("Flows", m_flows, true););
- }
-
- template
- bool network_flow::choose_leaving_edge() {
- TRACE("network_flow", tout << "choose_leaving_edge...\n";);
- node src = m_graph.get_source(m_enter_id);
- node tgt = m_graph.get_target(m_enter_id);
- m_delta.set_invalid();
- edge_id leave_id;
- svector path;
- svector against;
- m_tree->get_path(src, tgt, path, against);
- SASSERT(path.size() >= 1);
- for (unsigned i = 0; i < path.size(); ++i) {
- edge_id e_id = path[i];
- if (against[i] && (!m_delta || m_flows[e_id] < *m_delta)) {
- m_delta = m_flows[e_id];
- leave_id = e_id;
- }
- }
-
- if (m_delta) {
- m_leave_id = leave_id;
-
- TRACE("network_flow", {
- tout << "Found leaving edge " << m_leave_id;
- tout << " between node " << m_graph.get_source(m_leave_id);
- tout << " and node " << m_graph.get_target(m_leave_id) << " with delta = " << *m_delta << "...\n";
- });
- return true;
- }
- TRACE("network_flow", tout << "Can't find a leaving edge... The problem is unbounded.\n";);
- return false;
- }
-
- template
- void network_flow::update_spanning_tree() {
- m_tree->update(m_enter_id, m_leave_id);
- }
-
- template
- bool network_flow::choose_entering_edge(pivot_rule pr) {
- pivot_rule_impl * pivot;
- switch (pr) {
- case FIRST_ELIGIBLE:
- pivot = alloc(first_eligible_pivot, m_graph, m_potentials, m_states, m_enter_id);
- break;
- case BEST_ELIGIBLE:
- pivot = alloc(best_eligible_pivot, m_graph, m_potentials, m_states, m_enter_id);
- break;
- case CANDIDATE_LIST:
- pivot = alloc(candidate_list_pivot, m_graph, m_potentials, m_states, m_enter_id);
- break;
- default:
- UNREACHABLE();
- }
- return pivot->choose_entering_edge();
- }
-
- // Minimize cost flows
- template
- min_flow_result network_flow::min_cost(pivot_rule pr) {
- initialize();
- while (choose_entering_edge(pr)) {
- bool bounded = choose_leaving_edge();
- if (!bounded) return UNBOUNDED;
- update_flows();
- if (m_enter_id != m_leave_id) {
- SASSERT(edge_in_tree(m_leave_id));
- SASSERT(!edge_in_tree(m_enter_id));
- m_states[m_enter_id] = BASIS;
- m_states[m_leave_id] = LOWER;
- update_spanning_tree();
- update_potentials();
- TRACE("network_flow", tout << "Spanning tree:\n"; display_spanning_tree(tout););
- SASSERT(check_well_formed());
- }
- }
- if (is_infeasible()) return INFEASIBLE;
- TRACE("network_flow", tout << "Found optimal solution.\n";);
- SASSERT(check_optimal());
- return OPTIMAL;
- }
-
- template
- bool network_flow::is_infeasible() {
-#if 0
- // Flows of artificial arcs should be zero
- unsigned num_nodes = m_graph.get_num_nodes();
- unsigned num_edges = m_graph.get_num_edges();
- for (unsigned i = 1; i < num_nodes; ++i) {
- if (m_flows[num_edges - i].is_pos()) return true;
- }
-#endif
- return false;
- }
-
- // Get the optimal solution
- template
- typename network_flow::numeral network_flow::get_optimal_solution(vector & result, bool is_dual) {
- numeral objective_value = numeral::zero();
- unsigned num_edges = m_graph.get_num_edges();
- for (unsigned i = 0; i < num_edges; ++i) {
- if (m_states[i] == BASIS)
- {
- objective_value += m_flows[i].get_rational() * m_graph.get_weight(i);
- }
- }
- result.reset();
- if (is_dual) {
- result.append(m_potentials);
- }
- else {
- result.append(m_flows);
- }
- return objective_value;
- }
-
- template
- bool network_flow::edge_in_tree(edge_id id) const {
- return m_states[id] == BASIS;
- }
-
- template
- bool network_flow::check_well_formed() {
- SASSERT(m_tree->check_well_formed());
- SASSERT(!m_delta || !(*m_delta).is_neg());
-
- // m_flows are zero on non-basic edges
- for (unsigned i = 0; i < m_flows.size(); ++i) {
- SASSERT(!m_flows[i].is_neg());
- SASSERT(m_states[i] == BASIS || m_flows[i].is_zero());
- }
-
- unsigned num_edges = m_graph.get_num_edges();
- for (unsigned i = 0; i < num_edges; ++i) {
- if (m_states[i] == BASIS) {
- dl_var src = m_graph.get_source(i);
- dl_var tgt = m_graph.get_target(i);
- numeral weight = m_graph.get_weight(i);
- SASSERT(m_potentials[src] - m_potentials[tgt] == weight);
- }
- }
-
- return true;
- }
-
- template
- bool network_flow::check_optimal() {
- numeral total_cost = numeral::zero();
- unsigned num_edges = m_graph.get_num_edges();
- for (unsigned i = 0; i < num_edges; ++i) {
- if (m_states[i] == BASIS) {
- total_cost += m_flows[i].get_rational() * m_graph.get_weight(i);
- }
- }
-
- for (unsigned i = 0; i < num_edges; ++i) {
- dl_var src = m_graph.get_source(i);
- dl_var tgt = m_graph.get_target(i);
- numeral weight = m_graph.get_weight(i);
- SASSERT(m_potentials[src] - m_potentials[tgt] <= weight);
- }
-
- // m_flows are zero on non-basic edges
- for (unsigned i = 0; i < m_flows.size(); ++i) {
- SASSERT(m_states[i] == BASIS || m_flows[i].is_zero());
- }
- numeral total_balance = numeral::zero();
- for (unsigned i = 0; i < m_potentials.size(); ++i) {
- total_balance += m_balances[i] * m_potentials[i];
- }
- std::cout << "Total balance: " << total_balance << ", total cost: " << total_cost << std::endl;
- return total_cost == total_balance;
- }
-
template
void network_flow::display_spanning_tree(std::ofstream & os) {
- ++m_step;;
- std::string prefix = "T";
- prefix.append(std::to_string(m_step));
- prefix.append("_");
- unsigned root = m_graph.get_num_nodes() - 1;
- for (unsigned i = 0; i < root; ++i) {
- os << prefix << i << "[shape=circle,label=\"" << prefix << i << " [";
- os << m_potentials[i] << "/" << m_balances[i] << "]\"];\n";
- }
- os << prefix << root << "[shape=doublecircle,label=\"" << prefix << root << " [";
- os << m_potentials[root] << "/" << m_balances[root] << "]\"];\n";
+ ++m_step;;
+ std::string prefix = "T";
+ prefix.append(std::to_string(m_step));
+ prefix.append("_");
+ unsigned root = m_graph.get_num_nodes() - 1;
+ for (unsigned i = 0; i < root; ++i) {
+ os << prefix << i << "[shape=circle,label=\"" << prefix << i << " [";
+ os << m_potentials[i] << "/" << m_balances[i] << "]\"];\n";
+ }
+ os << prefix << root << "[shape=doublecircle,label=\"" << prefix << root << " [";
+ os << m_potentials[root] << "/" << m_balances[root] << "]\"];\n";
- unsigned num_edges = m_graph.get_num_edges();
- for (unsigned i = 0; i < num_edges; ++i) {
- os << prefix << m_graph.get_source(i) << " -> " << prefix << m_graph.get_target(i);
- if (m_states[i] == BASIS) {
- os << "[color=red,penwidth=3.0,label=\"" << m_flows[i] << "/" << m_graph.get_weight(i) << "\"];\n";
- }
- else {
- os << "[label=\"" << m_flows[i] << "/" << m_graph.get_weight(i) << "\"];\n";
- }
+ unsigned num_edges = m_graph.get_num_edges();
+ for (unsigned i = 0; i < num_edges; ++i) {
+ os << prefix << m_graph.get_source(i) << " -> " << prefix << m_graph.get_target(i);
+ if (edge_in_tree(i)) {
+ os << "[color=red,penwidth=3.0,label=\"" << m_flows[i] << "/" << m_graph.get_weight(i) << "\"];\n";
}
- os << std::endl;
+ else {
+ os << "[label=\"" << m_flows[i] << "/" << m_graph.get_weight(i) << "\"];\n";
+ }
+ }
+ os << std::endl;
}
}
diff --git a/src/smt/spanning_tree_base.h b/src/smt/spanning_tree_base.h
index f48ce8f4f..3d3ee7b3a 100644
--- a/src/smt/spanning_tree_base.h
+++ b/src/smt/spanning_tree_base.h
@@ -25,15 +25,8 @@ Notes:
namespace smt {
template
- inline std::string pp_vector(std::string const & label, TV v, bool has_header = false) {
+ inline std::string pp_vector(std::string const & label, TV v) {
std::ostringstream oss;
- if (has_header) {
- oss << "Index ";
- for (unsigned i = 0; i < v.size(); ++i) {
- oss << i << " ";
- }
- oss << std::endl;
- }
oss << label << " ";
for (unsigned i = 0; i < v.size(); ++i) {
oss << v[i] << " ";
diff --git a/src/smt/spanning_tree_def.h b/src/smt/spanning_tree_def.h
index 60bbde291..5dbcd2b8b 100644
--- a/src/smt/spanning_tree_def.h
+++ b/src/smt/spanning_tree_def.h
@@ -50,7 +50,7 @@ namespace smt {
}
TRACE("network_flow", {
- tout << pp_vector("Predecessors", m_pred, true) << pp_vector("Threads", m_thread);
+ tout << pp_vector("Predecessors", m_pred) << pp_vector("Threads", m_thread);
tout << pp_vector("Depths", m_depth) << pp_vector("Tree", m_tree);
});
}
@@ -204,7 +204,7 @@ namespace smt {
SASSERT(in_subtree_t2(v));
TRACE("network_flow", {
- tout << pp_vector("Predecessors", m_pred, true) << pp_vector("Threads", m_thread);
+ tout << pp_vector("Predecessors", m_pred) << pp_vector("Threads", m_thread);
tout << pp_vector("Depths", m_depth) << pp_vector("Tree", m_tree);
});
}