diff --git a/src/smt/diff_logic.h b/src/smt/diff_logic.h index 2717e4a92..74a5dbbd2 100644 --- a/src/smt/diff_logic.h +++ b/src/smt/diff_logic.h @@ -931,6 +931,24 @@ public: return found; } + // Return true if there is an edge source --> target. + // If there is such edge, return it in parameter e. + bool get_edge(dl_var source, dl_var target, edge & e) { + edge_id_vector & edges = m_out_edges[source]; + typename edge_id_vector::iterator it = edges.begin(); + typename edge_id_vector::iterator end = edges.end(); + bool found = false; + for (; it != end; ++it) { + edge_id e_id = *it; + edge & e0 = m_edges[e_id]; + if (e0.is_enabled() && e0.get_target() == target && !found) { + e = e0; + found = true; + } + } + return found; + } + template void enumerate_edges(dl_var source, dl_var target, Functor& f) { edge_id_vector & edges = m_out_edges[source]; diff --git a/src/smt/network_flow.h b/src/smt/network_flow.h new file mode 100644 index 000000000..d10ccf005 --- /dev/null +++ b/src/smt/network_flow.h @@ -0,0 +1,102 @@ +/*++ +Copyright (c) 2013 Microsoft Corporation + +Module Name: + + network_flow.h + +Abstract: + + Implements Network Simplex algorithm for min cost flow problem + +Author: + + Anh-Dung Phan (t-anphan) 2013-10-24 + +Notes: + + This will be used to solve the dual of min cost flow problem + i.e. optimization of difference constraint. + + We need a function to reduce DL constraints to min cost flow problem + and another function to convert from min cost flow solution to DL solution. + + It remains unclear how to convert DL assignment to a basic feasible solution of Network Simplex. + A naive approach is to run an algorithm on max flow in order to get a spanning tree. + + The network_simplex class hasn't had multiple pivoting strategies yet. + +--*/ +#ifndef _NETWORK_FLOW_H_ +#define _NETWORK_FLOW_H_ + +#include"inf_rational.h" +#include"diff_logic.h" + +namespace smt { + + // Solve minimum cost flow problem using Network Simplex algorithm + template + class network_flow : private Ext { + struct GExt : public Ext { + typedef literal explanation; + }; + + typedef dl_var node; + typedef dl_edge edge; + typedef dl_graph graph; + typedef typename Ext::numeral numeral; + graph m_graph; + + // Denote supply/demand b_i on node i + vector m_balances; + + vector m_potentials; + + // Keep optimal solution of the min cost flow problem + inf_rational m_objective; + + // Basic feasible flows + vector m_flows; + + // Denote the spanning tree of basic edges + vector m_basics; + // Denote non-basic edges with flow 0 for uncapicitated networks + vector m_nonbasics; + + // Store the parent of a node in the spanning tree + svector m_pred; + // Store the number of edge on the path to the root + svector m_depth; + // Store the pointer to the next node in depth first search ordering + svector m_thread; + + public: + // Initialize the network with a feasible spanning tree + void initialize(); + + void compute_potentials(); + + void compute_flows(); + + // If all reduced costs are non-negative, the current flow is optimal + // If not optimal, return a violating edge in the corresponding variable + bool is_optimal(edge & violating_edge); + + // Send as much flow as possible around the cycle, the first basic edge with flow 0 will leave + edge choose_leaving_edge(const edge & entering_edge); + + void update_basics(const edge & entering_edge, const edge & leaving_edge); + + bool is_unbounded(); + + // Compute the optimal solution + void get_optimal_solution(numeral & objective, vector & flows); + + // Minimize cost flows + // Return true if found an optimal solution, and return false if unbounded + bool min_cost(); + }; +} + +#endif diff --git a/src/smt/network_flow_def.h b/src/smt/network_flow_def.h new file mode 100644 index 000000000..8ce0a2458 --- /dev/null +++ b/src/smt/network_flow_def.h @@ -0,0 +1,148 @@ +/*++ +Copyright (c) 2013 Microsoft Corporation + +Module Name: + + network_flow_def.h + +Abstract: + + Implements Network Simplex algorithm for min cost flow problem + +Author: + + Anh-Dung Phan (t-anphan) 2013-10-24 + +Notes: + +--*/ + +#ifndef _NETWORK_FLOW_DEF_H_ +#define _NETWORK_FLOW_DEF_H_ + +#include"network_flow.h" + +namespace smt { + + template + void network_flow::initialize() { + // TODO: construct an initial spanning tree i.e. inializing m_pred, m_depth and m_thread. + compute_potentials(); + compute_flows(); + } + + template + void network_flow::compute_potentials() { + SASSERT(!m_potentials.empty()); + SASSERT(!m_thread.empty()); + SASSERT(m_thread.size() == m_pred.size()); + + array potentials; + std::copy(m_potentials.begin(), m_potentials.end(), potentials); + rational zero(0); + potentials[0] = zero; + node next = m_thread[0]; + + while (next != 0) { + node current = m_pred[next]; + edge e; + if (m_graph.get_edge(current, next, e)) { + potentials[next] = potentials[current] - e.get_weight(); + } + if (m_graph.get_edge(next, current, e)) { + potentials[next] = potentials[current] + e.get_weight(); + } + next = m_thread[next]; + } + std::copy(potentials.begin(), potentials.end(), m_potentials); + } + + template + void network_flow::compute_flows() { + vector balances(m_balances); + numeral zero; + m_flows.fill(zero); + vector basics(m_basics); + // TODO: need a way to find a leaf node of a spanning tree + while (!basics.empty()) { + return; + } + } + + template + bool network_flow::is_optimal(edge & violating_edge) { + typename vector::iterator it = m_nonbasics.begin(); + typename vector::iterator end = m_nonbasics.end(); + bool found = false; + for (unsigned int i = 0; i < m_nonbasics.size(); ++i) { + edge & e = m_nonbasics[i]; + if (e.is_enabled()) { + node source = e.get_source(); + node target = e.get_target(); + numeral cost = e.get_weight() - m_potentials[source] + m_potentials[target]; + // Choose the first negative-cost edge to be the violating edge + // TODO: add multiple pivoting strategies + if (cost < 0) { + violating_edge = e; + found = true; + break; + } + } + } + return !found; + } + + template + dl_edge::GExt> network_flow::choose_leaving_edge(const edge & entering_edge) { + node source = entering_edge.get_source(); + node target = entering_edge.get_target(); + while (source != target) { + if (m_depth[source] > m_depth[target]) + source = m_pred[source]; + else if (m_depth[source] < m_depth[target]) + target = m_pred[target]; + else { + source = m_pred[source]; + target = m_pred[target]; + } + } + edge e; + m_graph.get_edge(source, target, e); + return e; + } + + template + void network_flow::update_basics(const edge & entering_edge, const edge & leaving_edge) { + + } + + template + bool network_flow::is_unbounded() { + return false; + } + + // Get the optimal solution + template + void network_flow::get_optimal_solution(numeral & objective, vector & flows) { + flows.reset(); + flows.append(m_flows); + // TODO: calculate objective value + } + + // Minimize cost flows + // Return true if found an optimal solution, and return false if unbounded + template + bool network_flow::min_cost() { + initialize(); + edge & entering_edge; + while (!is_optimal(entering_edge)) { + edge & leaving_edge = choose_leaving_edge(); + update_tree(entering_edge, leaving_edge); + if (is_unbounded()) + return false; + } + return true; + } +} + +#endif diff --git a/src/smt/theory_diff_logic.h b/src/smt/theory_diff_logic.h index f904a600a..9447b0d3e 100644 --- a/src/smt/theory_diff_logic.h +++ b/src/smt/theory_diff_logic.h @@ -38,6 +38,7 @@ Revision History: #include"numeral_factory.h" #include"smt_clause.h" #include"theory_opt.h" +#include"network_flow_def.h" // The DL theory can represent term such as n + k, where n is an enode and k is a numeral. namespace smt { @@ -312,6 +313,8 @@ namespace smt { void internalize_objective(app * n, objective_term & objective); + network_flow m_network_flow; + private: virtual void new_eq_eh(theory_var v1, theory_var v2, justification& j); @@ -390,75 +393,6 @@ namespace smt { srdl_ext() : m_epsilon(s_integer(0),true) {} }; - // Solve minimum cost flow problem using Network Simplex algorithm - class network_simplex { - svector m_nodes; - svector m_edges; - // Denote costs c_ij on edge (i, j) - vector m_costs; - // Denote supply/demand b_i on node i - vector m_potentials; - - // Keep optimal solution of the min cost flow problem - inf_rational m_objective; - - // Data structure of spanning trees - svector m_pred; - svector m_depth; - svector m_thread; - - public: - // Initialize the network with a feasible spanning tree - virtual void initialize(); - - virtual void compute_potentials(); - - virtual void compute_flows(); - - // If all reduced costs are non-negative, the current flow is optimal - // If not optimal, return a violated edge in the corresponding variable - virtual bool is_optimal(edge_id & violated_edge); - - // Send as much flow as possible around the cycle, the first basic edge with flow 0 will leave - virtual edge_id choose_leaving_edge(); - - virtual void update_tree(edge_id entering_edge, edge_id leaving_edge); - - virtual bool is_unbounded(); - - // Compute the optimal solution - virtual void compute_solution(); - - // Minimize cost flows - // Return true if found an optimal solution, and return false if unbounded - bool minimize() { - initialize(); - compute_potentials(); - compute_flows(); - edge_id entering_edge; - while (!is_optimal(entering_edge)) { - edge_id leaving_edge = choose_leaving_edge(); - update_tree(entering_edge, leaving_edge); - - if (is_unbounded()) - return false; - } - return true; - } - }; - - /* Notes: - - We need a function to reduce DL constraints to min cost flow problem - and another function to convert from min cost flow solution to DL solution. - - It remains unclear how to convert DL assignment to a basic feasible solution of Network Simplex. - A naive approach is to run an algorithm on max flow in order to get a spanning tree. - - The network_simplex class hasn't had multiple pivoting strategies yet. - */ - - typedef theory_diff_logic theory_idl; typedef theory_diff_logic theory_fidl; typedef theory_diff_logic theory_rdl;