mirror of
https://github.com/Z3Prover/z3
synced 2025-04-22 16:45:31 +00:00
Update Network Simplex implementation
This commit is contained in:
Anh-Dung Phan
parent
d78d22deb6
commit
906bbb4eeb
4 changed files with 200 additions and 152 deletions
|
|
@ -598,7 +598,7 @@ namespace datalog {
|
|||
// 0 <= y - x - k - 1
|
||||
if (is_le(to_app(cond->get_arg(0)), x, k, y, is_int) && is_int) {
|
||||
k.neg();
|
||||
k -= rational::one();
|
||||
k -= rational::one();
|
||||
std::swap(x, y);
|
||||
return true;
|
||||
}
|
||||
|
|
|
|||
|
|
@ -38,6 +38,10 @@ namespace smt {
|
|||
// Solve minimum cost flow problem using Network Simplex algorithm
|
||||
template<typename Ext>
|
||||
class network_flow : private Ext {
|
||||
enum edge_state {
|
||||
NON_BASIS = 0,
|
||||
BASIS = 1
|
||||
};
|
||||
typedef dl_var node;
|
||||
typedef dl_edge<Ext> edge;
|
||||
typedef dl_graph<Ext> graph;
|
||||
|
|
@ -46,56 +50,66 @@ namespace smt {
|
|||
graph m_graph;
|
||||
|
||||
// Denote supply/demand b_i on node i
|
||||
vector<numeral> m_balances;
|
||||
vector<fin_numeral> m_balances;
|
||||
|
||||
// Duals of flows which are convenient to compute dual solutions
|
||||
vector<numeral> m_potentials;
|
||||
|
||||
// Keep optimal solution of the min cost flow problem
|
||||
inf_int_rational m_objective;
|
||||
numeral m_objective_value;
|
||||
|
||||
// Costs on edges
|
||||
vector<fin_numeral> const & m_costs;
|
||||
vector<fin_numeral> m_costs;
|
||||
|
||||
// Basic feasible flows
|
||||
vector<numeral> m_flows;
|
||||
|
||||
svector<edge_state> m_states;
|
||||
|
||||
// An element is true if the corresponding edge points upwards (compared to the root node)
|
||||
svector<bool> m_upwards;
|
||||
|
||||
// Store the parent of a node in the spanning tree
|
||||
// Store the parent of a node i in the spanning tree
|
||||
svector<node> m_pred;
|
||||
// Store the number of edge on the path to the root
|
||||
// Store the number of edge on the path from node i to the root
|
||||
svector<int> m_depth;
|
||||
// Store the pointer to the next node in depth first search ordering
|
||||
// Store the pointer from node i to the next node in depth first search ordering
|
||||
svector<node> m_thread;
|
||||
// Reverse orders of m_thread
|
||||
svector<node> m_rev_thread;
|
||||
// Store a final node of the sub tree rooted at node i
|
||||
svector<node> m_final;
|
||||
// Number of nodes in the sub tree rooted at node i
|
||||
svector<int> m_num_node;
|
||||
|
||||
bool m_is_optimal;
|
||||
edge_id m_entering_edge;
|
||||
edge_id m_leaving_edge;
|
||||
node m_join_node;
|
||||
numeral m_delta;
|
||||
|
||||
public:
|
||||
|
||||
network_flow(graph & g, vector<fin_numeral> const & costs);
|
||||
network_flow(graph & g, vector<fin_numeral> const & balances);
|
||||
|
||||
// Initialize the network with a feasible spanning tree
|
||||
void initialize();
|
||||
|
||||
void compute_potentials();
|
||||
void update_potentials();
|
||||
|
||||
void compute_flows();
|
||||
void update_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_id & violating_edge);
|
||||
// If all reduced costs are non-negative, return false since the current spanning tree is optimal
|
||||
// Otherwise return true and update m_entering_edge
|
||||
bool choose_entering_edge();
|
||||
|
||||
// Send as much flow as possible around the cycle, the first basic edge with flow 0 will leave
|
||||
edge_id choose_leaving_edge(edge_id entering_edge);
|
||||
// Return false if the problem is unbounded
|
||||
bool choose_leaving_edge();
|
||||
|
||||
void update_spanning_tree(edge_id entering_edge, edge_id leaving_edge);
|
||||
void update_spanning_tree();
|
||||
|
||||
bool is_unbounded();
|
||||
|
||||
// Compute the optimal solution
|
||||
void get_optimal_solution(numeral & objective, vector<numeral> & flows);
|
||||
numeral get_optimal_solution(vector<numeral> & result, bool is_dual);
|
||||
|
||||
// Minimize cost flows
|
||||
// Return true if found an optimal solution, and return false if unbounded
|
||||
|
|
|
|||
|
|
@ -25,170 +25,212 @@ Notes:
|
|||
namespace smt {
|
||||
|
||||
template<typename Ext>
|
||||
network_flow<Ext>::network_flow(graph & g, vector<fin_numeral> const& costs) :
|
||||
network_flow<Ext>::network_flow(graph & g, vector<fin_numeral> const& balances) :
|
||||
m_graph(g),
|
||||
m_costs(costs) {
|
||||
m_balances(balances) {
|
||||
unsigned num_nodes = m_balances.size() + 1;
|
||||
unsigned num_edges = m_graph.get_num_edges();
|
||||
|
||||
vector<edge> const & es = m_graph.get_all_edges();
|
||||
for (unsigned i = 0; i < num_edges; ++i) {
|
||||
fin_numeral cost(es[i].get_weight().get_rational());
|
||||
m_costs.push_back(cost);
|
||||
}
|
||||
|
||||
m_balances.resize(num_nodes);
|
||||
for (unsigned i = 0; i < balances.size(); ++i) {
|
||||
m_costs.push_back(balances[i]);
|
||||
}
|
||||
|
||||
m_potentials.resize(num_nodes);
|
||||
m_costs.resize(num_edges);
|
||||
m_flows.resize(num_edges);
|
||||
m_states.resize(num_edges);
|
||||
|
||||
m_upwards.resize(num_nodes);
|
||||
m_pred.resize(num_nodes);
|
||||
m_depth.resize(num_nodes);
|
||||
m_thread.resize(num_nodes);
|
||||
m_rev_thread.resize(num_nodes);
|
||||
m_final.resize(num_nodes);
|
||||
m_num_node.resize(num_nodes);
|
||||
}
|
||||
|
||||
template<typename Ext>
|
||||
void network_flow<Ext>::initialize() {
|
||||
// TODO: construct an initial spanning tree i.e. inializing m_pred, m_depth and m_thread.
|
||||
compute_potentials();
|
||||
compute_flows();
|
||||
}
|
||||
// Create an artificial root node to construct initial spanning tree
|
||||
unsigned num_nodes = m_balances.size();
|
||||
unsigned num_edges = m_graph.get_num_edges();
|
||||
node root = num_nodes;
|
||||
m_pred[root] = -1;
|
||||
m_thread[root] = 0;
|
||||
m_rev_thread[0] = root;
|
||||
m_num_node[root] = num_nodes + 1;
|
||||
m_final[root] = root - 1;
|
||||
m_potentials[root] = numeral::zero();
|
||||
|
||||
template<typename Ext>
|
||||
void network_flow<Ext>::compute_potentials() {
|
||||
SASSERT(!m_potentials.empty());
|
||||
SASSERT(!m_thread.empty());
|
||||
SASSERT(m_thread.size() == m_pred.size());
|
||||
|
||||
m_potentials.set(0, numeral::zero());
|
||||
node target = m_thread[0];
|
||||
|
||||
while (target != 0) {
|
||||
node source = m_pred[target];
|
||||
edge_id e_id;
|
||||
if (m_graph.get_edge_id(source, target, e_id)) {
|
||||
m_potentials.set(target, m_potentials[source] - m_graph.get_weight(e_id));
|
||||
}
|
||||
if (m_graph.get_edge_id(target, source, e_id)) {
|
||||
m_potentials.set(target, m_potentials[source] + m_graph.get_weight(e_id));
|
||||
}
|
||||
target = m_thread[target];
|
||||
fin_numeral sum_supply = fin_numeral::zero();
|
||||
for (unsigned i = 0; i < m_balances.size(); ++i) {
|
||||
sum_supply += m_balances[i];
|
||||
}
|
||||
m_balances[root] = -sum_supply;
|
||||
|
||||
m_states.resize(num_nodes + num_edges);
|
||||
m_states.fill(NON_BASIS);
|
||||
|
||||
// Create artificial edges and initialize the spanning tree
|
||||
for (unsigned i = 0; i < num_nodes; ++i) {
|
||||
m_upwards[i] = m_balances[i] >= fin_numeral::zero();
|
||||
m_pred[i] = root;
|
||||
m_depth[i] = 1;
|
||||
m_thread[i] = i + 1;
|
||||
m_final[i] = i;
|
||||
m_rev_thread[i] = (i = 0) ? root : i - 1;
|
||||
m_num_node[i] = 1;
|
||||
m_states[num_edges + i] = BASIS;
|
||||
}
|
||||
}
|
||||
|
||||
template<typename Ext>
|
||||
void network_flow<Ext>::compute_flows() {
|
||||
vector<numeral> balances(m_balances);
|
||||
|
||||
// OPTIMIZE: Need a set data structure for efficiently removing elements
|
||||
vector<edge_id> basics;
|
||||
while (!basics.empty()) {
|
||||
// Find a leaf node of a spanning tree
|
||||
node target;
|
||||
for (unsigned int i = 0; i < m_thread.size(); ++i) {
|
||||
if (m_depth[i] <= m_depth[m_thread[i]]) {
|
||||
target = i;
|
||||
break;
|
||||
}
|
||||
}
|
||||
node source = m_pred[target];
|
||||
void network_flow<Ext>::update_potentials() {
|
||||
node src = m_graph.get_source(m_entering_edge);
|
||||
node tgt = m_graph.get_source(m_entering_edge);
|
||||
numeral cost = m_graph.get_weight(m_entering_edge);
|
||||
numeral change = m_upwards[src] ? (cost - m_potentials[src] + m_potentials[tgt]) :
|
||||
(-cost + m_potentials[src] - m_potentials[tgt]);
|
||||
node last = m_thread[m_final[src]];
|
||||
for (node u = src; u != last; u = m_thread[u]) {
|
||||
m_potentials[u] += change;
|
||||
}
|
||||
}
|
||||
|
||||
template<typename Ext>
|
||||
void network_flow<Ext>::update_flows() {
|
||||
numeral val = m_state[m_entering_edge] == NON_BASIS ? numeral::zero() : m_delta;
|
||||
m_flows[m_entering_edge] += val;
|
||||
for (unsigned u = m_graph.get_source(m_entering_edge); u != m_join_node; u = m_pred[u]) {
|
||||
edge_id e_id;
|
||||
if (m_graph.get_edge_id(source, target, e_id)) {
|
||||
m_flows.set(e_id, -m_graph.get_weight(basics[target]));
|
||||
basics[source] += basics[target];
|
||||
basics.erase(e_id);
|
||||
}
|
||||
else if (m_graph.get_edge_id(target, source, e_id)) {
|
||||
m_flows.set(e_id, m_graph.get_weight(basics[target]));
|
||||
basics[source] += basics[target];
|
||||
basics.erase(e_id);
|
||||
}
|
||||
m_graph.get_edge_id(u, m_pred[u], e_id);
|
||||
m_flows[e_id] += m_upwards[u] ? -val : val;
|
||||
}
|
||||
for (unsigned u = m_graph.get_target(m_entering_edge); u != m_join_node; u = m_pred[u]) {
|
||||
edge_id e_id;
|
||||
m_graph.get_edge_id(u, m_pred[u], e_id);
|
||||
m_flows[e_id] += m_upwards[u] ? val : -val;
|
||||
}
|
||||
}
|
||||
|
||||
template<typename Ext>
|
||||
bool network_flow<Ext>::is_optimal(edge_id & violating_edge) {
|
||||
// TODO: how to get nonbasics vector?
|
||||
vector<edge> nonbasics;
|
||||
typename vector<edge>::iterator it = nonbasics.begin();
|
||||
typename vector<edge>::iterator end = nonbasics.end();
|
||||
bool found = false;
|
||||
for (unsigned int i = 0; i < nonbasics.size(); ++i) {
|
||||
edge & e = nonbasics[i];
|
||||
if (e.is_enabled()) {
|
||||
bool network_flow<Ext>::choose_entering_edge() {
|
||||
vector<edge> const & es = m_graph.get_all_edges();
|
||||
for (unsigned int i = 0; i < es.size(); ++i) {
|
||||
edge const & e = es[i];
|
||||
edge_id e_id;
|
||||
if (e.is_enabled() && m_graph.get_edge_id(e.get_source(), e.get_target(), e_id) && m_states[e_id] == BASIS) {
|
||||
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
|
||||
numeral zero(0);
|
||||
if (cost < zero) {
|
||||
edge_id e_id;
|
||||
m_graph.get_edge_id(source, target, e_id);
|
||||
violating_edge = e_id;
|
||||
found = true;
|
||||
break;
|
||||
if (cost < numeral::zero()) {
|
||||
m_entering_edge = e_id;
|
||||
return true;
|
||||
}
|
||||
}
|
||||
}
|
||||
return !found;
|
||||
return false;
|
||||
}
|
||||
|
||||
template<typename Ext>
|
||||
edge_id network_flow<Ext>::choose_leaving_edge(edge_id entering_edge) {
|
||||
node source = m_graph.get_source(entering_edge);
|
||||
node target = m_graph.get_target(entering_edge);
|
||||
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];
|
||||
bool network_flow<Ext>::choose_leaving_edge() {
|
||||
node source = m_graph.get_source(m_entering_edge);
|
||||
node target = m_graph.get_target(m_entering_edge);
|
||||
node u = source, v = target;
|
||||
while (u != v) {
|
||||
if (m_depth[u] > m_depth[v])
|
||||
u = m_pred[u];
|
||||
else if (m_depth[u] < m_depth[v])
|
||||
v = m_pred[v];
|
||||
else {
|
||||
source = m_pred[source];
|
||||
target = m_pred[target];
|
||||
u = m_pred[u];
|
||||
v = m_pred[v];
|
||||
}
|
||||
}
|
||||
edge_id e_id;
|
||||
m_graph.get_edge_id(source, target, e_id);
|
||||
return e_id;
|
||||
}
|
||||
|
||||
template<typename Ext>
|
||||
void network_flow<Ext>::update_spanning_tree(edge_id entering_edge, edge_id leaving_edge) {
|
||||
// Need special handling in case two edges are identical
|
||||
SASSERT(entering_edge != leaving_edge);
|
||||
|
||||
// Update potentials
|
||||
node target = m_upwards[leaving_edge] ? m_graph.get_source(leaving_edge) : m_graph.get_target(leaving_edge);
|
||||
numeral src_pot = m_potentials[m_graph.get_source(entering_edge)];
|
||||
numeral tgt_pot = m_potentials[m_graph.get_target(entering_edge)];
|
||||
numeral weight = m_graph.get_weight(entering_edge);
|
||||
numeral change = m_upwards[entering_edge] ? (weight - src_pot + tgt_pot) : (-weight + src_pot - tgt_pot);
|
||||
m_potentials[target] += change;
|
||||
node start = m_thread[target];
|
||||
while (m_depth[start] > m_depth[target]) {
|
||||
m_potentials[start] += change;
|
||||
start = m_thread[start];
|
||||
// Found first common ancestor of source and target
|
||||
m_join_node = u;
|
||||
// FIXME: need to get truly finite value
|
||||
numeral infty = numeral(INT_MAX);
|
||||
m_delta = infty;
|
||||
node src, tgt;
|
||||
// Send flows along the path from source to the ancestor
|
||||
for (unsigned u = source; u != m_join_node; u = m_pred[u]) {
|
||||
edge_id e_id;
|
||||
m_graph.get_edge_id(u, m_pred[u], e_id);
|
||||
numeral d = m_upwards[u] ? m_flows[e_id] : infty;
|
||||
if (d < m_delta) {
|
||||
m_delta = d;
|
||||
src = u;
|
||||
tgt = m_pred[u];
|
||||
}
|
||||
}
|
||||
|
||||
// Send flows along the path from target to the ancestor
|
||||
for (unsigned u = target; u != m_join_node; u = m_pred[u]) {
|
||||
edge_id e_id;
|
||||
m_graph.get_edge_id(u, m_pred[u], e_id);
|
||||
numeral d = m_upwards[u] ? infty : m_flows[e_id];
|
||||
if (d <= m_delta) {
|
||||
m_delta = d;
|
||||
src = u;
|
||||
tgt = m_pred[u];
|
||||
}
|
||||
}
|
||||
|
||||
if (m_delta < infty) {
|
||||
m_graph.get_edge_id(src, tgt, m_leaving_edge);
|
||||
return true;
|
||||
}
|
||||
return false;
|
||||
}
|
||||
|
||||
template<typename Ext>
|
||||
bool network_flow<Ext>::is_unbounded() {
|
||||
return false;
|
||||
void network_flow<Ext>::update_spanning_tree() {
|
||||
|
||||
}
|
||||
|
||||
// Get the optimal solution
|
||||
template<typename Ext>
|
||||
void network_flow<Ext>::get_optimal_solution(numeral & objective, vector<numeral> & flows) {
|
||||
SASSERT(m_is_optimal);
|
||||
flows.reset();
|
||||
flows.append(m_flows);
|
||||
objective = numeral::zero();
|
||||
for (unsigned int i = 0; i < m_flows.size(); ++i) {
|
||||
objective += m_costs[i] * m_flows[i];
|
||||
typename network_flow<Ext>::numeral network_flow<Ext>::get_optimal_solution(vector<numeral> & result, bool is_dual) {
|
||||
m_objective_value = numeral::zero();
|
||||
for (unsigned i = 0; i < m_flows.size(); ++i) {
|
||||
m_objective_value += m_costs[i] * m_flows[i];
|
||||
}
|
||||
result.reset();
|
||||
if (is_dual) {
|
||||
result.append(m_potentials);
|
||||
}
|
||||
else {
|
||||
result.append(m_flows);
|
||||
}
|
||||
return m_objective_value;
|
||||
}
|
||||
|
||||
// Minimize cost flows
|
||||
// Return true if found an optimal solution, and return false if unbounded
|
||||
template<typename Ext>
|
||||
bool network_flow<Ext>::min_cost() {
|
||||
SASSERT(!m_graph.get_all_edges().empty());
|
||||
initialize();
|
||||
edge_id entering_edge;
|
||||
while (!is_optimal(entering_edge)) {
|
||||
edge_id leaving_edge = choose_leaving_edge(entering_edge);
|
||||
update_spanning_tree(entering_edge, leaving_edge);
|
||||
if (is_unbounded()) {
|
||||
m_is_optimal = false;
|
||||
return m_is_optimal;
|
||||
while (choose_entering_edge()) {
|
||||
bool bounded = choose_leaving_edge();
|
||||
if (!bounded) return false;
|
||||
if (m_entering_edge != m_leaving_edge) {
|
||||
m_states[m_entering_edge] = BASIS;
|
||||
m_states[m_leaving_edge] = NON_BASIS;
|
||||
update_spanning_tree();
|
||||
update_potentials();
|
||||
}
|
||||
}
|
||||
m_is_optimal = true;
|
||||
return m_is_optimal;
|
||||
return true;
|
||||
}
|
||||
}
|
||||
|
||||
|
|
|
|||
|
|
@ -1016,34 +1016,26 @@ bool theory_diff_logic<Ext>::maximize(theory_var v) {
|
|||
// m_graph as well.
|
||||
dl_graph<GExt> g;
|
||||
vector<dl_edge<GExt> > const& es = m_graph.get_all_edges();
|
||||
dl_var offset = m_graph.get_num_edges();
|
||||
for (unsigned i = 0; i < es.size(); ++i) {
|
||||
dl_edge<GExt> const & e = es[i];
|
||||
if (e.is_enabled()) {
|
||||
g.enable_edge(g.add_edge(e.get_source(), e.get_target(), e.get_weight(), e.get_explanation()));
|
||||
g.enable_edge(g.add_edge(e.get_target() + offset, e.get_source() + offset, e.get_weight(), e.get_explanation()));
|
||||
}
|
||||
}
|
||||
|
||||
// Objective coefficients now become costs
|
||||
vector<fin_numeral> base_costs, aux_costs;
|
||||
// Objective coefficients now become balances
|
||||
vector<fin_numeral> balances;
|
||||
for (unsigned i = 0; i < objective.size(); ++i) {
|
||||
fin_numeral cost(objective[i].second);
|
||||
base_costs.push_back(cost);
|
||||
aux_costs.push_back(-cost);
|
||||
fin_numeral balance(objective[i].second);
|
||||
balances.push_back(balance);
|
||||
}
|
||||
vector<fin_numeral> costs;
|
||||
costs.append(base_costs);
|
||||
costs.append(aux_costs);
|
||||
|
||||
network_flow<GExt> net_flow(g, costs);
|
||||
network_flow<GExt> net_flow(g, balances);
|
||||
bool is_optimal = net_flow.min_cost();
|
||||
if (is_optimal) {
|
||||
numeral objective_value;
|
||||
vector<numeral> flows;
|
||||
net_flow.get_optimal_solution(objective_value, flows);
|
||||
m_objective_value = objective_value.get_rational();
|
||||
// TODO: return the model of the optimal solution
|
||||
vector<numeral> potentials;
|
||||
m_objective_value = net_flow.get_optimal_solution(potentials, true);
|
||||
// TODO: return the model of the optimal solution from potential
|
||||
}
|
||||
return is_optimal;
|
||||
}
|
||||
|
|
|
|||
Loading…
Add table
Add a link
Reference in a new issue