diff --git a/src/smt/diff_logic.h b/src/smt/diff_logic.h
index 74a5dbbd2..71b638b9e 100644
--- a/src/smt/diff_logic.h
+++ b/src/smt/diff_logic.h
@@ -932,21 +932,23 @@ public:
     }
 
     // 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) {
+    // If there is such edge, return its edge_id in parameter id.
+    bool get_edge_id(dl_var source, dl_var target, edge_id & id) {
         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;
+            id = *it;
+            edge & e = m_edges[id];
+            if (e.is_enabled() && e.get_target() == target) {
+                return true;
             }
         }
-        return found;
+        return false;
+    }
+
+    edges const & get_all_edges() const {
+        return m_edges;
     }
 
     template<typename Functor>
diff --git a/src/smt/network_flow.h b/src/smt/network_flow.h
index d10ccf005..81f9aff1f 100644
--- a/src/smt/network_flow.h
+++ b/src/smt/network_flow.h
@@ -38,31 +38,30 @@ 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;        
+        typedef dl_edge<Ext> edge;
+        typedef dl_graph<Ext> graph;     
+        typedef typename Ext::numeral numeral;
+        typedef typename Ext::fin_numeral fin_numeral;
         graph m_graph;
 
         // Denote supply/demand b_i on node i
         vector<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_rational m_objective;
+        inf_int_rational m_objective;
+
+        // Costs on edges
+        vector<fin_numeral> const & m_costs;
 
         // 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;
+        // 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
         svector<node> m_pred;
@@ -71,7 +70,12 @@ namespace smt {
         // Store the pointer to the next node in depth first search ordering
         svector<node> m_thread;
 
+        bool m_is_optimal;
+
     public:
+
+        network_flow(graph & g, vector<fin_numeral> const & costs);
+
         // Initialize the network with a feasible spanning tree
         void initialize();
 
@@ -81,13 +85,13 @@ namespace smt {
             
         // 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);
+        bool is_optimal(edge_id & 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);
+        edge_id choose_leaving_edge(edge_id entering_edge);
 
+        void update_spanning_tree(edge_id entering_edge, edge_id leaving_edge);
+        
         bool is_unbounded();
 
         // Compute the optimal solution
diff --git a/src/smt/network_flow_def.h b/src/smt/network_flow_def.h
index 8ce0a2458..464747141 100644
--- a/src/smt/network_flow_def.h
+++ b/src/smt/network_flow_def.h
@@ -24,6 +24,12 @@ Notes:
 
 namespace smt {
 
+    template<typename Ext>
+    network_flow<Ext>::network_flow(graph & g, vector<fin_numeral> const& costs) :
+        m_graph(g),
+        m_costs(costs) {
+    }
+
     template<typename Ext>
     void network_flow<Ext>::initialize() {
         // TODO: construct an initial spanning tree i.e. inializing m_pred, m_depth and m_thread.
@@ -36,54 +42,73 @@ namespace smt {
         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();
+        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(next, current, e)) {
-                potentials[next] = potentials[current] + e.get_weight();
+            if (m_graph.get_edge_id(target, source, e_id)) {
+                m_potentials.set(target, m_potentials[source] + m_graph.get_weight(e_id));
             }
-            next = m_thread[next];
+            target = m_thread[target];
         }
-        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
+        
+        // OPTIMIZE: Need a set data structure for efficiently removing elements
+        vector<edge_id> basics;
         while (!basics.empty()) {
-            return;
+            // 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];
+            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);
+            }
         }
     }
     
     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 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 < m_nonbasics.size(); ++i) {
-            edge & e = m_nonbasics[i];
+        for (unsigned int i = 0; i < nonbasics.size(); ++i) {
+            edge & e = 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;
+                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;
                 }
@@ -93,9 +118,9 @@ namespace smt {
     }
 
     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();
+    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];
@@ -106,14 +131,28 @@ namespace smt {
                 target = m_pred[target];
             }
         }
-        edge e;
-        m_graph.get_edge(source, target, e);
-        return e;
+        edge_id e_id;
+        m_graph.get_edge_id(source, target, e_id);
+        return e_id;
     }
 
     template<typename Ext>
-    void network_flow<Ext>::update_basics(const edge & entering_edge, const edge & leaving_edge) {
-
+    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];
+        }
     }
 
     template<typename Ext>
@@ -124,24 +163,34 @@ namespace smt {
     // 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);
-        // TODO: calculate objective value
+        numeral cost(0);
+        for (unsigned int i = 0; i < m_flows.size(); ++i) {
+            // FIXME: this * operator is not supported
+            cost += m_costs[i] * m_flows[i];
+        }
+        objective = cost;
     }
 
     // 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 & entering_edge;
+        edge_id entering_edge;
         while (!is_optimal(entering_edge)) {
-            edge & leaving_edge = choose_leaving_edge();
-            update_tree(entering_edge, leaving_edge);
-            if (is_unbounded())
-                return false;
+            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;
+            }
         }
-        return true;
+        m_is_optimal = true;
+        return m_is_optimal;
     }
 }
 
diff --git a/src/smt/theory_diff_logic.h b/src/smt/theory_diff_logic.h
index ab47b9423..f7f574cd0 100644
--- a/src/smt/theory_diff_logic.h
+++ b/src/smt/theory_diff_logic.h
@@ -38,7 +38,7 @@ Revision History:
 #include"numeral_factory.h"
 #include"smt_clause.h"
 #include"theory_opt.h"
-#include"network_flow_def.h"
+#include"network_flow.h"
 
 // The DL theory can represent term such as n + k, where n is an enode and k is a numeral.
 namespace smt {
@@ -307,14 +307,13 @@ namespace smt {
         virtual bool maximize(theory_var v);
         virtual theory_var add_objective(app* term);
         virtual inf_eps_rational<inf_rational> get_objective_value(theory_var v);
+        numeral m_objective_value;
 
         typedef vector <std::pair<theory_var, rational> > objective_term;
         vector<objective_term> m_objectives;
-        vector<rational>       m_objective_vars;
+        vector<rational>       m_objective_consts;
 
-        void internalize_objective(app * n, objective_term & objective);
-
-        network_flow<Ext> m_network_flow;
+        bool internalize_objective(expr * n, rational const& m, rational& r, objective_term & objective);
 
     private:        
 
@@ -368,6 +367,7 @@ namespace smt {
         // TODO: It doesn't need to be a rational, but a bignum integer.
         static const bool m_int_theory = true;
         typedef rational numeral;
+        typedef rational fin_numeral;
         numeral     m_epsilon;
         idl_ext() : m_epsilon(1) {}
     };
@@ -376,6 +376,7 @@ namespace smt {
         // TODO: It doesn't need to be a rational, but a bignum integer.
         static const bool m_int_theory = true;
         typedef s_integer numeral;
+        typedef s_integer fin_numeral;
         numeral m_epsilon;
         sidl_ext() : m_epsilon(1) {}
     };
@@ -383,13 +384,15 @@ namespace smt {
     struct rdl_ext {
         static const bool m_int_theory = false;
         typedef inf_int_rational numeral;
-        numeral      m_epsilon;
+        typedef rational fin_numeral;
+        numeral      m_epsilon;        
         rdl_ext() : m_epsilon(rational(), true) {}
     };
 
     struct srdl_ext {
         static const bool m_int_theory = false;
         typedef inf_s_integer numeral;
+        typedef s_integer fin_numeral;
         numeral m_epsilon;
         srdl_ext() : m_epsilon(s_integer(0),true) {}
     };
diff --git a/src/smt/theory_diff_logic_def.h b/src/smt/theory_diff_logic_def.h
index f392ac423..478b80964 100644
--- a/src/smt/theory_diff_logic_def.h
+++ b/src/smt/theory_diff_logic_def.h
@@ -28,6 +28,7 @@ Revision History:
 #include"ast_pp.h"
 #include"warning.h"
 #include"smt_model_generator.h"
+#include"network_flow_def.h"
 
 using namespace smt;
 
@@ -55,8 +56,6 @@ std::ostream& theory_diff_logic<Ext>::atom::display(theory_diff_logic const& th,
 template<typename Ext>
 void theory_diff_logic<Ext>::nc_functor::reset() {
     m_antecedents.reset();
-    m_objectives.reset();
-    m_objective_vars.reset();
 }
 
 
@@ -777,6 +776,8 @@ void theory_diff_logic<Ext>::reset_eh() {
     m_agility               = 0.5;
     m_is_lia                = true;
     m_non_diff_logic_exprs  = false;
+    m_objectives      .reset();
+    m_objective_consts.reset();
     theory::reset_eh();
 }
 
@@ -1000,14 +1001,51 @@ void theory_diff_logic<Ext>::get_implied_bound_antecedents(edge_id bridge_edge,
 
 template<typename Ext>
 bool theory_diff_logic<Ext>::maximize(theory_var v) {
+    objective_term const& objective = m_objectives[v];
+
     IF_VERBOSE(1,
-               objective_term const& objective = m_objectives[v];
                for (unsigned i = 0; i < objective.size(); ++i) {
                    verbose_stream() << "coefficient " << objective[i].second << " of theory_var " << objective[i].first << "\n";
                }
-               verbose_stream() << m_objective_vars[v] << "\n";);
+               verbose_stream() << m_objective_consts[v] << "\n";);
     NOT_IMPLEMENTED_YET();
-    return false;
+    // Double the number of edges in the new graph
+
+    // NSB review: what about disabled edges? They should not be added, right?
+    // For efficiency we probably want to reuse m_graph and keep extra edges on the side or add them to 
+    // 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;
+    for (unsigned i = 0; i < objective.size(); ++i) {
+        fin_numeral cost(objective[i].second);
+        base_costs.push_back(cost);
+        aux_costs.push_back(-cost);
+    }
+    vector<fin_numeral> costs;
+    costs.append(base_costs);
+    costs.append(aux_costs);
+    
+    network_flow<GExt> net_flow(g, costs);
+    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        
+    }
+    return is_optimal;
 }
 
 template<typename Ext>
@@ -1015,59 +1053,55 @@ theory_var theory_diff_logic<Ext>::add_objective(app* term) {
     objective_term objective;
     theory_var result = m_objectives.size();
     rational q(1), r(0);
-    if (!internalize_objective(term, q, r, objective)) {
-        result = null_theory_var;
+    if (internalize_objective(term, q, r, objective)) {
+        m_objectives.push_back(objective);
+        m_objective_consts.push_back(r);
     }
     else {
-        m_objectives.push_back(objective);
-        m_objective_vars.push_back(r);
+        result = null_theory_var;
     }
     return result; 
 }
 
 template<typename Ext>
 inf_eps_rational<inf_rational> theory_diff_logic<Ext>::get_objective_value(theory_var v) { 
+    NOT_IMPLEMENTED_YET();
     inf_rational objective;
     inf_eps_rational<inf_rational> val(objective);     
     return val; 
 }
 
 template<typename Ext>
-bool theory_diff_logic<Ext>::internalize_objective(app * n, rational& q, objective_term & objective) {
+bool theory_diff_logic<Ext>::internalize_objective(expr * n, rational const& m, rational& q, objective_term & objective) {
 
     // Compile term into objective_term format
     rational r;
+    expr* x, *y;
     if (m_util.is_numeral(n, r)) {
         q += r;
     }
     else if (m_util.is_add(n)) {
-        for (unsigned i = 0; i < n->get_num_args(); ++i) {
-            if (!internalize_objective(to_app(n->get_arg(i)), objective)) {
+        for (unsigned i = 0; i < to_app(n)->get_num_args(); ++i) {
+            if (!internalize_objective(to_app(n)->get_arg(i), m, q, objective)) {
                 return false;
             }
-        };
+        }
     }
-    else if (m_util.is_mul(n)) {
-        SASSERT(n->get_num_args() == 2);
-        rational r;
-        app * lhs = to_app(n->get_arg(0));
-        app * rhs = to_app(n->get_arg(1));
-        SASSERT(m_util.is_numeral(lhs, r) || m_util.is_numeral(rhs, r));
-
-        if (!m_util.is_numeral(lhs, r))
-            std::swap(lhs, rhs);
-        
-        m_util.is_numeral(lhs, r);
-        theory_var v = mk_var(rhs);
-        objective.push_back(std::make_pair(v, r));
+    else if (m_util.is_mul(n, x, y) && m_util.is_numeral(x, r)) {
+        return internalize_objective(y, m*r, q, objective);
     }
-    else if (n->get_family_id() == m_util.get_family_id()) {
+    else if (m_util.is_mul(n, y, x) && m_util.is_numeral(x, r)) {
+        return internalize_objective(y, m*r, q, objective);
+    }
+    else if (!is_app(n)) {
+        return false;
+    }
+    else if (to_app(n)->get_family_id() == m_util.get_family_id()) {
         return false;
     }
     else {
-        theory_var v = mk_var(n);
-        rational r(1);
-        objective.push_back(std::make_pair(v, r));
+        theory_var v = mk_var(to_app(n));
+        objective.push_back(std::make_pair(v, m));
     }
     return true;
 }
diff --git a/src/util/inf_int_rational.h b/src/util/inf_int_rational.h
index c9879f0ee..b1b1fb89f 100644
--- a/src/util/inf_int_rational.h
+++ b/src/util/inf_int_rational.h
@@ -155,6 +155,17 @@ class inf_int_rational {
 	return *this; 
     }
 
+    inf_int_rational & operator*=(const rational & r) { 
+        if (!r.is_int32()) {
+            throw default_exception("multiplication with large rational is not possible");
+        }
+        m_first  *= r;
+        m_second *= r.get_int32();
+	return *this; 
+    }
+
+
+
     inf_int_rational & operator-=(const inf_int_rational & r) { 
         m_first  -= r.m_first;
         m_second -= r.m_second;
@@ -344,6 +355,10 @@ inline inf_int_rational operator+(const inf_int_rational & r1, const inf_int_rat
     return inf_int_rational(r1) += r2; 
 }
 
+inline inf_int_rational operator*(const rational & r1, const inf_int_rational & r2) { 
+    return inf_int_rational(r2) *= r1;
+}
+
 inline inf_int_rational operator-(const inf_int_rational & r1, const inf_int_rational & r2) { 
     return inf_int_rational(r1) -= r2; 
 }
diff --git a/src/util/rational.h b/src/util/rational.h
index fc03228bf..cf3c7b238 100644
--- a/src/util/rational.h
+++ b/src/util/rational.h
@@ -111,6 +111,11 @@ public:
         return INT_MIN <= v && v <= INT_MAX;
     }
 
+    int get_int32() const {
+        SASSERT(is_int32());
+        return (int)get_int64();
+    }
+
     double get_double() const { return m().get_double(m_val); }
 
     rational const & get_rational() const { return *this; }