diff --git a/src/muz/rel/dl_interval_relation.cpp b/src/muz/rel/dl_interval_relation.cpp
index c04269b02..88c262cbe 100644
--- a/src/muz/rel/dl_interval_relation.cpp
+++ b/src/muz/rel/dl_interval_relation.cpp
@@ -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;
             }
diff --git a/src/smt/network_flow.h b/src/smt/network_flow.h
index 81f9aff1f..37afda136 100644
--- a/src/smt/network_flow.h
+++ b/src/smt/network_flow.h
@@ -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
diff --git a/src/smt/network_flow_def.h b/src/smt/network_flow_def.h
index 7f6add1eb..f39042bac 100644
--- a/src/smt/network_flow_def.h
+++ b/src/smt/network_flow_def.h
@@ -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;
     }
 }
 
diff --git a/src/smt/theory_diff_logic_def.h b/src/smt/theory_diff_logic_def.h
index 478b80964..9fa97097f 100644
--- a/src/smt/theory_diff_logic_def.h
+++ b/src/smt/theory_diff_logic_def.h
@@ -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;
 }