diff --git a/src/opt/maxsmt.cpp b/src/opt/maxsmt.cpp
index d27d17790..cc1b1ba42 100644
--- a/src/opt/maxsmt.cpp
+++ b/src/opt/maxsmt.cpp
@@ -105,10 +105,15 @@ namespace opt {
         return r;
     }
 
-    void maxsmt::update_lower(rational const& r) {
-        if (m_lower > r)  m_lower = r;
+    void maxsmt::update_lower(rational const& r, bool override) {
+        if (m_lower > r || override)  m_lower = r;
     }
 
+    void maxsmt::update_upper(rational const& r, bool override) {
+        if (m_upper < r || override)  m_upper = r;
+    }
+    
+
     void maxsmt::get_model(model_ref& mdl) {
         mdl = m_model.get();
     }
diff --git a/src/opt/maxsmt.h b/src/opt/maxsmt.h
index c4cf33cac..cd26c40e5 100644
--- a/src/opt/maxsmt.h
+++ b/src/opt/maxsmt.h
@@ -73,7 +73,8 @@ namespace opt {
         rational get_value() const;
         rational get_lower() const;
         rational get_upper() const;
-        void update_lower(rational const& r);
+        void update_lower(rational const& r, bool override);
+        void update_upper(rational const& r, bool override);
         void get_model(model_ref& mdl);
         bool get_assignment(unsigned index) const;
         void display_answer(std::ostream& out) const;        
diff --git a/src/opt/opt_context.cpp b/src/opt/opt_context.cpp
index a16654614..4f2a5eba2 100644
--- a/src/opt/opt_context.cpp
+++ b/src/opt/opt_context.cpp
@@ -152,7 +152,7 @@ namespace opt {
         IF_VERBOSE(1, verbose_stream() << "(optimize:sat)\n";);
         s.get_model(m_model);
         m_optsmt.setup(s);
-        update_lower();
+        update_lower(true);
         switch (m_objectives.size()) {
         case 0:
             return is_sat;
@@ -234,91 +234,147 @@ namespace opt {
         return r;
     }
 
+    class context::pareto : public pareto_callback {
+        context& ctx;
+        ast_manager& m;
+        expr_ref mk_ge(expr* t, expr* s) {
+            expr_ref result(m);
+            if (ctx.m_bv.is_bv(t)) {
+                result = ctx.m_bv.mk_ule(s, t);
+            }
+            else {
+                result = ctx.m_arith.mk_ge(t, s);
+            }
+            return result;
+        }
+    public:
+        pareto(context& ctx):ctx(ctx),m(ctx.m) {}
+
+
+        virtual void yield(model_ref& mdl) {
+            ctx.m_model = mdl;
+            ctx.update_lower(true);
+            for (unsigned i = 0; i < ctx.m_objectives.size(); ++i) {
+                objective const& obj = ctx.m_objectives[i];
+                switch(obj.m_type) {
+                case O_MINIMIZE:
+                case O_MAXIMIZE:
+                    ctx.m_optsmt.update_upper(obj.m_index, ctx.m_optsmt.get_lower(obj.m_index), true);
+                    break;
+                case O_MAXSMT: {
+                    rational r = ctx.m_maxsmts.find(obj.m_id)->get_lower();
+                    ctx.m_maxsmts.find(obj.m_id)->update_upper(r, true);
+                    break;
+                }
+                }
+            }
+
+            IF_VERBOSE(1, ctx.display_assignment(verbose_stream()););
+        }
+        virtual unsigned num_objectives() {
+            return ctx.m_objectives.size();
+        }
+        virtual expr_ref mk_le(unsigned i, model_ref& mdl) {
+            objective const& obj = ctx.m_objectives[i];
+            expr_ref val(m), result(m), term(m);
+            mk_term_val(mdl, obj, term, val);
+            switch (obj.m_type) {
+            case O_MINIMIZE:
+                result = mk_ge(term, val);
+                break;
+            case O_MAXSMT:
+                result = mk_ge(term, val);
+                break;
+            case O_MAXIMIZE:
+                result = mk_ge(val, term);
+                break;
+            }
+            return result;
+        }
+        virtual expr_ref mk_ge(unsigned i, model_ref& mdl) {
+            objective const& obj = ctx.m_objectives[i];
+            expr_ref val(m), result(m), term(m);
+            mk_term_val(mdl, obj, term, val);
+            switch (obj.m_type) {
+            case O_MINIMIZE:
+                result = mk_ge(val, term);
+                break;
+            case O_MAXSMT:
+                result = mk_ge(val, term);
+                break;
+            case O_MAXIMIZE:
+                result = mk_ge(term, val);
+                break;
+            }
+            return result;
+        }
+
+        virtual expr_ref mk_gt(unsigned i, model_ref& mdl) {
+            expr_ref result = mk_le(i, mdl);
+            result = m.mk_not(result);
+            return result;
+        }
+    private:
+        void mk_term_val(model_ref& mdl, objective const& obj, expr_ref& term, expr_ref& val) {
+            rational r;
+            switch (obj.m_type) {
+            case O_MINIMIZE:
+            case O_MAXIMIZE:
+                term = obj.m_term;
+                break;
+            case O_MAXSMT: {
+                unsigned sz = obj.m_terms.size();
+                expr_ref_vector sum(m);
+                expr_ref zero(m);
+                zero = ctx.m_arith.mk_numeral(rational(0), false);
+                for (unsigned i = 0; i < sz; ++i) {
+                    expr* t = obj.m_terms[i];
+                    rational const& w = obj.m_weights[i];
+                    sum.push_back(m.mk_ite(t, ctx.m_arith.mk_numeral(w, false), zero));                    
+                }
+                if (sum.empty()) {
+                    term = zero;
+                }
+                else {
+                    term = ctx.m_arith.mk_add(sum.size(), sum.c_ptr());
+                }                           
+                break;
+            }
+            }
+            VERIFY(mdl->eval(term, val) && ctx.is_numeral(val, r));
+        }        
+
+#if 0
+        // use PB
+        expr_ref mk_pb_cmp(objective const& obj, model_ref& mdl, bool is_ge) {
+            rational r, sum(0);
+            expr_ref val(m), result(m);
+            unsigned sz = obj.m_terms.size();
+            pb_util pb(m);
+
+            for (unsigned i = 0; i < sz; ++i) {
+                expr* t = obj.m_terms[i];
+                VERIFY(mdl->eval(t, val));
+                if (m.is_true(val)) {
+                    sum += r;
+                }
+            }
+            if (is_ge) {
+                result = pb.mk_ge(obj.m_terms.size(), obj.m_weights.c_ptr(), obj.m_terms.c_ptr(), r);
+            }
+            else {
+                result = pb.mk_le(obj.m_terms.size(), obj.m_weights.c_ptr(), obj.m_terms.c_ptr(), r);
+            }
+        }
+#endif
+    };
 
     lbool context::execute_pareto() {
-        opt_solver& s = get_solver();        
-        expr_ref val(m);
-        rational r;
-        lbool is_sat = l_true;
-        vector<bounds_t> bounds;
-        for (unsigned i = 0; i < m_objectives.size(); ++i) {
-            objective const& obj = m_objectives[i];
-            if (obj.m_type == O_MAXSMT) {
-                IF_VERBOSE(0, verbose_stream() << "Pareto optimization is not supported for MAXSMT\n";);
-                return l_undef;
-            }
-            solver::scoped_push _sp(s);
-            is_sat = m_optsmt.pareto(obj.m_index);
-            if (is_sat != l_true) {
-                return is_sat;
-            }
-            if (!m_optsmt.get_upper(obj.m_index).is_finite()) {
-                return l_undef;
-            }
-            bounds_t bound;            
-            for (unsigned j = 0; j < m_objectives.size(); ++j) {
-                objective const& obj_j = m_objectives[j];
-                inf_eps lo = m_optsmt.get_lower(obj_j.m_index);
-                inf_eps hi = m_optsmt.get_upper(obj_j.m_index);
-                bound.push_back(std::make_pair(lo, hi));
-            }            
-            bounds.push_back(bound);
-        }
-        for (unsigned i = 0; i < bounds.size(); ++i) {
-            for (unsigned j = 0; j < bounds.size(); ++j) {
-                objective const& obj = m_objectives[j];
-                bounds[i][j].second = bounds[j][j].second;                
-            }
-            IF_VERBOSE(0, display_bounds(verbose_stream() << "new bound\n", bounds[i]););
-        }
-
-        for (unsigned i = 0; i < bounds.size(); ++i) {
-            bounds_t b = bounds[i];
-            vector<inf_eps> mids;
-            solver::scoped_push _sp(s);
-            for (unsigned j = 0; j < b.size(); ++j) {
-                objective const& obj = m_objectives[j];
-                inf_eps mid = (b[j].second - b[j].first)/rational(2);
-                mids.push_back(mid);
-                expr_ref ge = s.mk_ge(obj.m_index, mid);            
-                s.assert_expr(ge);
-            }
-            is_sat = execute_box();
-            switch(is_sat) {
-            case l_undef: 
-                return is_sat;
-            case l_true: {
-                bool at_bound = true; 
-                for (unsigned j = 0; j < b.size(); ++j) {
-                    objective const& obj = m_objectives[j];
-                    if (m_model->eval(obj.m_term, val) && is_numeral(val, r)) {
-                        mids[j] = inf_eps(r);
-                    }
-                    at_bound = at_bound && mids[j] == b[j].second;
-                    b[j].second = mids[j];
-                }
-                IF_VERBOSE(0, display_bounds(verbose_stream() << "new bound\n", b););
-                if (!at_bound) {
-                    bounds.push_back(b);
-                }
-                break;
-            }
-            case l_false: {
-                bounds_t b2(b);
-                for (unsigned j = 0; j < b.size(); ++j) {
-                    b2[j].second = mids[j];
-                    if (j > 0) {
-                        b2[j-1].second = b[j-1].second;
-                    }
-                    IF_VERBOSE(0, display_bounds(verbose_stream() << "refined bound\n", b2););
-                    bounds.push_back(b2);
-                }
-                break;
-            }
-            }
-        }
-        
-        return is_sat;
+        pareto cb(*this);
+        m_pareto = alloc(gia_pareto, m, cb, m_solver.get(), m_params);
+        return (*(m_pareto.get()))();
+        // NB. stack reference cb is out of scope after return.
+        // NB. fix race condition for set_cancel
     }
 
     void context::display_bounds(std::ostream& out, bounds_t const& b) const {
@@ -652,7 +708,7 @@ namespace opt {
         }
     }
 
-    void context::update_lower() {
+    void context::update_lower(bool override) {
         expr_ref val(m);
         rational r(0);
         for (unsigned i = 0; i < m_objectives.size(); ++i) {
@@ -660,12 +716,12 @@ namespace opt {
             switch(obj.m_type) {
             case O_MINIMIZE:
                 if (m_model->eval(obj.m_term, val) && is_numeral(val, r)) {
-                    m_optsmt.update_lower(obj.m_index, -r);
+                    m_optsmt.update_lower(obj.m_index, inf_eps(-r), override);
                 }
                 break;
             case O_MAXIMIZE:
                 if (m_model->eval(obj.m_term, val) && is_numeral(val, r)) {
-                    m_optsmt.update_lower(obj.m_index, r);
+                    m_optsmt.update_lower(obj.m_index, inf_eps(r), override);
                 }
                 break;
             case O_MAXSMT: {
@@ -681,7 +737,7 @@ namespace opt {
                     }
                 }
                 if (ok) {
-                    m_maxsmts.find(obj.m_id)->update_lower(r);
+                    m_maxsmts.find(obj.m_id)->update_lower(r, override);
                 }
                 break;
             }
@@ -816,6 +872,9 @@ namespace opt {
         if (m_simplify) {
             m_simplify->set_cancel(f);
         }
+        if (m_pareto) {
+            m_pareto->set_cancel(f);
+        }
         m_optsmt.set_cancel(f);
         map_t::iterator it = m_maxsmts.begin(), end = m_maxsmts.end();
         for (; it != end; ++it) {
diff --git a/src/opt/opt_context.h b/src/opt/opt_context.h
index 3fad9a59d..1912dac0d 100644
--- a/src/opt/opt_context.h
+++ b/src/opt/opt_context.h
@@ -14,19 +14,13 @@ Author:
 
 Notes:
 
-    TODO:
-
-    - type check objective term and assertions. It should pass basic sanity being
-      integer, real (, bit-vector) or other supported objective function type.
-
-    - add appropriate statistics tracking to opt::context
-
 --*/
 #ifndef _OPT_CONTEXT_H_
 #define _OPT_CONTEXT_H_
 
 #include "ast.h"
 #include "opt_solver.h"
+#include "opt_pareto.h"
 #include "optsmt.h"
 #include "maxsmt.h"
 #include "model_converter.h"
@@ -86,6 +80,7 @@ namespace opt {
         bv_util             m_bv;
         expr_ref_vector     m_hard_constraints;
         ref<opt_solver>     m_solver;
+        scoped_ptr<pareto_base>          m_pareto;
         params_ref          m_params;
         optsmt              m_optsmt;        
         map_t               m_maxsmts;
@@ -158,7 +153,7 @@ namespace opt {
         void from_fmls(expr_ref_vector const& fmls);
         void simplify_fmls(expr_ref_vector& fmls);
 
-        void update_lower();
+        void update_lower(bool override);
 
         inf_eps get_lower_as_num(unsigned idx);
         inf_eps get_upper_as_num(unsigned idx);
@@ -174,6 +169,8 @@ namespace opt {
 
         void validate_lex();
 
+        class pareto;
+
     };
 
 }
diff --git a/src/opt/opt_pareto.cpp b/src/opt/opt_pareto.cpp
new file mode 100644
index 000000000..6bffd19ef
--- /dev/null
+++ b/src/opt/opt_pareto.cpp
@@ -0,0 +1,197 @@
+/*++
+Copyright (c) 2014 Microsoft Corporation
+
+Module Name:
+
+    opt_pareto.cpp
+
+Abstract:
+
+    Pareto front utilities
+
+Author:
+
+    Nikolaj Bjorner (nbjorner) 2014-4-24
+
+Notes:
+
+   
+--*/
+
+#include "opt_pareto.h"
+#include "ast_pp.h"
+
+namespace opt {
+
+    // ---------------------
+    // GIA pareto algorithm
+   
+    lbool gia_pareto::operator()() {
+        model_ref model;
+        expr_ref fml(m);
+        lbool is_sat = m_solver->check_sat(0, 0);
+        while (is_sat == l_true) {
+            {
+                solver::scoped_push _s(*m_solver.get());
+                while (is_sat == l_true) {
+                    if (m_cancel) {
+                        return l_undef;
+                    }
+                    m_solver->get_model(model);
+                    // TBD: we can also use local search to tune solution coordinate-wise.
+                    mk_dominates(model);
+                    is_sat = m_solver->check_sat(0, 0);
+                }
+                if (is_sat == l_undef) {
+                    return l_undef;
+                }
+                is_sat = l_true;
+            }
+            cb.yield(model);
+            mk_not_dominated_by(model);
+            is_sat = m_solver->check_sat(0, 0);
+        }
+        if (is_sat == l_undef) {
+            return l_undef;
+        }
+        return l_true;
+    }
+
+    void pareto_base::mk_dominates(model_ref& model) {
+        unsigned sz = cb.num_objectives();
+        expr_ref fml(m);
+        expr_ref_vector gt(m), fmls(m);
+        for (unsigned i = 0; i < sz; ++i) {
+            fmls.push_back(cb.mk_ge(i, model));
+            gt.push_back(cb.mk_gt(i, model));
+        }
+        fmls.push_back(m.mk_or(gt.size(), gt.c_ptr()));
+        fml = m.mk_and(fmls.size(), fmls.c_ptr());
+        IF_VERBOSE(10, verbose_stream() << "dominates: " << fml << "\n";);
+        m_solver->assert_expr(fml);        
+    }
+
+    void pareto_base::mk_not_dominated_by(model_ref& model) {
+        unsigned sz = cb.num_objectives();
+        expr_ref fml(m);
+        expr_ref_vector le(m);
+        for (unsigned i = 0; i < sz; ++i) {
+            le.push_back(cb.mk_le(i, model));
+        }
+        fml = m.mk_not(m.mk_and(le.size(), le.c_ptr()));
+        IF_VERBOSE(10, verbose_stream() << "not dominated by: " << fml << "\n";);
+        m_solver->assert_expr(fml);        
+    }
+
+    // ---------------------------------
+    // OIA algorithm (without filtering)
+
+    lbool oia_pareto::operator()() {
+        model_ref model;
+        solver::scoped_push _s(*m_solver.get());
+        lbool is_sat = m_solver->check_sat(0, 0);
+        if (is_sat != l_true) {
+            return is_sat;
+        }
+        while (is_sat == l_true) {
+            if (m_cancel) {
+                return l_undef;
+            }
+            m_solver->get_model(model);
+            cb.yield(model);
+            mk_not_dominated_by(model);
+            is_sat = m_solver->check_sat(0, 0);
+        }
+        if (m_cancel) {
+            return l_undef;
+        }
+        return l_true;
+    }
+
+}
+
+#if 0
+        opt_solver& s = get_solver();        
+        expr_ref val(m);
+        rational r;
+        lbool is_sat = l_true;
+        vector<bounds_t> bounds;
+        for (unsigned i = 0; i < m_objectives.size(); ++i) {
+            objective const& obj = m_objectives[i];
+            if (obj.m_type == O_MAXSMT) {
+                IF_VERBOSE(0, verbose_stream() << "Pareto optimization is not supported for MAXSMT\n";);
+                return l_undef;
+            }
+            solver::scoped_push _sp(s);
+            is_sat = m_optsmt.pareto(obj.m_index);
+            if (is_sat != l_true) {
+                return is_sat;
+            }
+            if (!m_optsmt.get_upper(obj.m_index).is_finite()) {
+                return l_undef;
+            }
+            bounds_t bound;            
+            for (unsigned j = 0; j < m_objectives.size(); ++j) {
+                objective const& obj_j = m_objectives[j];
+                inf_eps lo = m_optsmt.get_lower(obj_j.m_index);
+                inf_eps hi = m_optsmt.get_upper(obj_j.m_index);
+                bound.push_back(std::make_pair(lo, hi));
+            }            
+            bounds.push_back(bound);
+        }
+        for (unsigned i = 0; i < bounds.size(); ++i) {
+            for (unsigned j = 0; j < bounds.size(); ++j) {
+                objective const& obj = m_objectives[j];
+                bounds[i][j].second = bounds[j][j].second;                
+            }
+            IF_VERBOSE(0, display_bounds(verbose_stream() << "new bound\n", bounds[i]););
+        }
+
+        for (unsigned i = 0; i < bounds.size(); ++i) {
+            bounds_t b = bounds[i];
+            vector<inf_eps> mids;
+            solver::scoped_push _sp(s);
+            for (unsigned j = 0; j < b.size(); ++j) {
+                objective const& obj = m_objectives[j];
+                inf_eps mid = (b[j].second - b[j].first)/rational(2);
+                mids.push_back(mid);
+                expr_ref ge = s.mk_ge(obj.m_index, mid);            
+                s.assert_expr(ge);
+            }
+            is_sat = execute_box();
+            switch(is_sat) {
+            case l_undef: 
+                return is_sat;
+            case l_true: {
+                bool at_bound = true; 
+                for (unsigned j = 0; j < b.size(); ++j) {
+                    objective const& obj = m_objectives[j];
+                    if (m_model->eval(obj.m_term, val) && is_numeral(val, r)) {
+                        mids[j] = inf_eps(r);
+                    }
+                    at_bound = at_bound && mids[j] == b[j].second;
+                    b[j].second = mids[j];
+                }
+                IF_VERBOSE(0, display_bounds(verbose_stream() << "new bound\n", b););
+                if (!at_bound) {
+                    bounds.push_back(b);
+                }
+                break;
+            }
+            case l_false: {
+                bounds_t b2(b);
+                for (unsigned j = 0; j < b.size(); ++j) {
+                    b2[j].second = mids[j];
+                    if (j > 0) {
+                        b2[j-1].second = b[j-1].second;
+                    }
+                    IF_VERBOSE(0, display_bounds(verbose_stream() << "refined bound\n", b2););
+                    bounds.push_back(b2);
+                }
+                break;
+            }
+            }
+        }
+        
+        return is_sat;
+#endif
diff --git a/src/opt/opt_pareto.h b/src/opt/opt_pareto.h
new file mode 100644
index 000000000..fa0243807
--- /dev/null
+++ b/src/opt/opt_pareto.h
@@ -0,0 +1,109 @@
+/*++
+Copyright (c) 2014 Microsoft Corporation
+
+Module Name:
+
+    opt_pareto.h
+
+Abstract:
+
+    Pareto front utilities
+
+Author:
+
+    Nikolaj Bjorner (nbjorner) 2014-4-24
+
+Notes:
+
+   
+--*/
+#ifndef _OPT_PARETO_H_
+#define _OPT_PARETO_H_
+
+#include "solver.h"
+#include "model.h"
+
+namespace opt {
+   
+    class pareto_callback {
+    public:
+        virtual void yield(model_ref& model) = 0;
+        virtual unsigned num_objectives() = 0;
+        virtual expr_ref mk_gt(unsigned i, model_ref& model) = 0;
+        virtual expr_ref mk_ge(unsigned i, model_ref& model) = 0;
+        virtual expr_ref mk_le(unsigned i, model_ref& model) = 0;
+    };
+    class pareto_base {
+    protected:
+        ast_manager&     m;
+        pareto_callback& cb;
+        volatile bool    m_cancel;
+        ref<solver>      m_solver;
+        params_ref       m_params;
+    public:
+        pareto_base(
+            ast_manager & m, 
+            pareto_callback& cb, 
+            solver* s, 
+            params_ref & p):
+            m(m),
+            cb(cb),            
+            m_cancel(false),
+            m_solver(s),
+            m_params(p) {
+        }
+        virtual ~pareto_base() {}
+        virtual void updt_params(params_ref & p) {
+            m_solver->updt_params(p);
+            m_params.copy(p);
+        }
+        virtual void collect_param_descrs(param_descrs & r) {
+            m_solver->collect_param_descrs(r);
+        }
+        virtual void collect_statistics(statistics & st) const {
+            m_solver->collect_statistics(st);
+        }
+        virtual void set_cancel(bool f) {
+            m_solver->set_cancel(f);
+            m_cancel = f;
+        }
+        virtual void display(std::ostream & out) const {
+            m_solver->display(out);
+        }
+        virtual lbool operator()() = 0;
+
+    protected:
+
+        void mk_dominates(model_ref& model);
+
+        void mk_not_dominated_by(model_ref& model);            
+    };
+    class gia_pareto : public pareto_base {
+    public:
+        gia_pareto(ast_manager & m, 
+                   pareto_callback& cb, 
+                   solver* s, 
+                   params_ref & p):
+            pareto_base(m, cb, s, p) {            
+        }
+        virtual ~gia_pareto() {}
+
+        virtual lbool operator()();
+    };
+
+    // opportunistic improvement algorithm.
+    class oia_pareto : public pareto_base {
+    public:
+        oia_pareto(ast_manager & m, 
+                   pareto_callback& cb, 
+                   solver* s, 
+                   params_ref & p):
+            pareto_base(m, cb, s, p) {            
+        }
+        virtual ~oia_pareto() {}
+
+        virtual lbool operator()();
+    };
+}
+
+#endif
diff --git a/src/opt/optsmt.cpp b/src/opt/optsmt.cpp
index 14e133a70..e5370e82d 100644
--- a/src/opt/optsmt.cpp
+++ b/src/opt/optsmt.cpp
@@ -115,11 +115,15 @@ namespace opt {
         return l_true;        
     }
 
-    void optsmt::update_lower(unsigned idx, rational const& r) {
-        inf_eps v(r);
-        if (m_lower[idx] < v) {
+    void optsmt::update_lower(unsigned idx, inf_eps const& v, bool override) {
+        if (m_lower[idx] < v || override) {
             m_lower[idx] = v;            
-            if (m_s) m_s->get_model(m_model);
+        }
+    }
+
+    void optsmt::update_upper(unsigned idx, inf_eps const& v, bool override) {
+        if (m_upper[idx] > v || override) {
+            m_upper[idx] = v;            
         }
     }
 
diff --git a/src/opt/optsmt.h b/src/opt/optsmt.h
index 060ae9aaa..5f4a2988f 100644
--- a/src/opt/optsmt.h
+++ b/src/opt/optsmt.h
@@ -60,7 +60,8 @@ namespace opt {
         inf_eps get_upper(unsigned index) const;
         void    get_model(model_ref& mdl);
 
-        void update_lower(unsigned idx, rational const& r);
+        void update_lower(unsigned idx, inf_eps const& r, bool override);
+        void update_upper(unsigned idx, inf_eps const& r, bool override);
 
         void reset();