WIP on min cost flow problem

Remarks:
1. Follow the template structure of diff_logic.h
2. Try to reuse dl_graph<Ext> with some ready-to-use graph algorithms
3. Need to add 'explanation' to 'GExt' in order to instantiate
dl_graph<_>

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<typename Functor>
     void enumerate_edges(dl_var source, dl_var target, Functor& f) {
         edge_id_vector & edges = m_out_edges[source];

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<typename Ext>
+    class network_flow : private Ext {
+        struct GExt : public Ext {
+            typedef literal explanation;
+        };
+
+        typedef dl_var node;
+        typedef dl_edge<GExt> edge;
+        typedef dl_graph<GExt> graph;     
+        typedef typename Ext::numeral numeral;        
+        graph m_graph;
+
+        // Denote supply/demand b_i on node i
+        vector<numeral> m_balances;
+
+        vector<numeral> m_potentials;
+
+        // Keep optimal solution of the min cost flow problem
+        inf_rational m_objective;
+
+        // Basic feasible flows
+        vector<numeral> m_flows;
+
+        // Denote the spanning tree of basic edges
+        vector<edge> m_basics;
+        // Denote non-basic edges with flow 0 for uncapicitated networks
+        vector<edge> m_nonbasics;
+
+        // Store the parent of a node in the spanning tree
+        svector<node> m_pred;
+        // Store the number of edge on the path to the root
+        svector<int> m_depth;
+        // Store the pointer to the next node in depth first search ordering
+        svector<node> 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<numeral> & flows);
+
+        // Minimize cost flows
+        // Return true if found an optimal solution, and return false if unbounded
+        bool min_cost();
+    };
+}
+
+#endif

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<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();
+    }
+
+    template<typename Ext>
+    void network_flow<Ext>::compute_potentials() {
+        SASSERT(!m_potentials.empty());
+        SASSERT(!m_thread.empty());
+        SASSERT(m_thread.size() == m_pred.size());
+
+        array<rational, m_potentials.size()> 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<typename Ext>
+    void network_flow<Ext>::compute_flows() {
+        vector<numeral> balances(m_balances);
+        numeral zero;
+        m_flows.fill(zero);
+        vector<edge> basics(m_basics);
+        // TODO: need a way to find a leaf node of a spanning tree
+        while (!basics.empty()) {
+            return;
+        }
+    }
+    
+    template<typename Ext>
+    bool network_flow<Ext>::is_optimal(edge & violating_edge) {        
+        typename vector<edge>::iterator it  = m_nonbasics.begin();
+        typename vector<edge>::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<typename Ext>
+    dl_edge<typename network_flow<Ext>::GExt> network_flow<Ext>::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<typename Ext>
+    void network_flow<Ext>::update_basics(const edge & entering_edge, const edge & leaving_edge) {
+
+    }
+
+    template<typename Ext>
+    bool network_flow<Ext>::is_unbounded() {
+        return false;
+    }
+
+    // Get the optimal solution
+    template<typename Ext>
+    void network_flow<Ext>::get_optimal_solution(numeral & objective, vector<numeral> & 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<typename Ext>
+    bool network_flow<Ext>::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

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<Ext> 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<dl_var>  m_nodes;
-        svector<edge_id> m_edges;
-        // Denote costs c_ij on edge (i, j)
-        vector<rational> m_costs;
-        // Denote supply/demand b_i on node i
-        vector<rational> m_potentials;
-
-        // Keep optimal solution of the min cost flow problem
-        inf_rational m_objective;
-
-        // Data structure of spanning trees
-        svector<dl_var> m_pred;
-        svector<int> m_depth;
-        svector<dl_var> 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<idl_ext>  theory_idl;
     typedef theory_diff_logic<sidl_ext> theory_fidl;
     typedef theory_diff_logic<rdl_ext> theory_rdl;