diff --git a/src/math/simplex/model_based_opt.cpp b/src/math/simplex/model_based_opt.cpp
index 86b3cfa0c..2b175506f 100644
--- a/src/math/simplex/model_based_opt.cpp
+++ b/src/math/simplex/model_based_opt.cpp
@@ -35,11 +35,113 @@ std::ostream& operator<<(std::ostream& out, opt::ineq_type ie) {
 
 namespace opt {
     
+    /**
+     * Convert a row ax + coeffs + coeff = 0 into a definition for x
+     *    x  = -(coeffs + coeff)/a 
+     * For rows ax + coeffs + coeff < 0 convert into
+     *    x  = -(coeffs + coeff - a)/a 
+     */
+    model_based_opt::def::def(row const& r, unsigned x) {
+        rational c = r.get_coefficient(x);
+        bool is_pos = c.is_pos();
+        for (var const & v : r.m_vars) {
+            if (v.m_id != x) { 
+                m_vars.push_back(v); 
+                if (is_pos) m_vars.back().m_coeff.neg(); 
+            }
+        }        
+        m_coeff = r.m_coeff;
+        if (is_pos) m_coeff.neg();
+        if (r.m_type == opt::t_lt) m_coeff += abs(c);
+        m_div = abs(c);
+        normalize();
+        SASSERT(m_div.is_pos());
+    }
+
+    model_based_opt::def model_based_opt::def::operator+(def const& other) const {
+        def result;
+        vector<var> const& vs1 = m_vars;
+        vector<var> const& vs2 = other.m_vars;
+        vector<var> & vs = result.m_vars;
+        rational c1(1), c2(1);
+        if (m_div != other.m_div) {
+            c1 = other.m_div;
+            c2 = m_div;
+        }        
+        unsigned i = 0, j = 0;
+        while (i < vs1.size() || j < vs2.size()) {
+            unsigned v1 = UINT_MAX, v2 = UINT_MAX;
+            if (i < vs1.size()) v1 = vs1[i].m_id;
+            if (j < vs2.size()) v2 = vs2[j].m_id;
+            if (v1 == v2) {
+                vs.push_back(vs1[i]);
+                vs.back().m_coeff *= c1; 
+                vs.back().m_coeff += c2 * vs2[j].m_coeff; 
+                ++i; ++j;
+                if (vs.back().m_coeff.is_zero()) {
+                    vs.pop_back();
+                }
+            }
+            else if (v1 < v2) {
+                vs.push_back(vs1[i]);
+                vs.back().m_coeff *= c1;                 
+            }
+            else {
+                vs.push_back(vs2[j]);
+                vs.back().m_coeff *= c2;                 
+            }
+        }
+        result.m_div = c1*m_div;
+        result.m_coeff = (m_coeff*c1) + (other.m_coeff*c2);
+        result.normalize();
+        return result;
+    }
+
+    model_based_opt::def model_based_opt::def::operator/(rational const& r) const {
+        def result(*this);
+        result.m_div *= r;
+        result.normalize();
+        return result;
+    }
+
+    model_based_opt::def model_based_opt::def::operator*(rational const& n) const {
+        def result(*this);
+        for (var& v : result.m_vars) {
+            v.m_coeff *= n;
+        }
+        result.m_coeff *= n;
+        result.normalize();
+        return result;
+    }
+
+    model_based_opt::def model_based_opt::def::operator+(rational const& n) const {
+        def result(*this);
+        result.m_coeff += n * result.m_div;
+        result.normalize();
+        return result;
+    }
+
+    void model_based_opt::def::normalize() {
+        if (m_div.is_one()) return;
+        rational g(m_div);
+        g = gcd(g, m_coeff);
+        for (var const& v : m_vars) {
+            g = gcd(g, abs(v.m_coeff));
+            if (g.is_one()) break;
+        }
+        if (!g.is_one()) {
+            for (var& v : m_vars) {
+                v.m_coeff /= g;
+            }
+            m_coeff /= g;
+            m_div /= g;
+        }
+    }
+
 
     model_based_opt::model_based_opt() {
         m_rows.push_back(row());
     }
-
         
     bool model_based_opt::invariant() {
         for (unsigned i = 0; i < m_rows.size(); ++i) {
@@ -105,14 +207,14 @@ namespace opt {
             if (find_bound(x, bound_row_index, bound_coeff, coeff.is_pos())) {
                 SASSERT(!bound_coeff.is_zero());
                 TRACE("opt", display(tout << "update: " << v << " ", objective());
-                      for (unsigned i = 0; i < m_above.size(); ++i) {
-                          display(tout << "resolve: ", m_rows[m_above[i]]);
+                      for (unsigned above : m_above) {
+                          display(tout << "resolve: ", m_rows[above]);
                       });
-                for (unsigned i = 0; i < m_above.size(); ++i) {
-                    resolve(bound_row_index, bound_coeff, m_above[i], x);
+                for (unsigned above : m_above) {
+                    resolve(bound_row_index, bound_coeff, above, x);
                 }
-                for (unsigned i = 0; i < m_below.size(); ++i) {
-                    resolve(bound_row_index, bound_coeff, m_below[i], x);
+                for (unsigned below : m_below) {
+                    resolve(bound_row_index, bound_coeff, below, x);
                 }
                 // coeff*x + objective <= ub
                 // a2*x + t2 <= 0
@@ -151,8 +253,7 @@ namespace opt {
         rational old_val = m_var2value[x];
         m_var2value[x] = val;
         unsigned_vector const& row_ids = m_var2row_ids[x];
-        for (unsigned i = 0; i < row_ids.size(); ++i) {            
-            unsigned row_id = row_ids[i];
+        for (unsigned row_id : row_ids) {
             rational coeff = get_coefficient(row_id, x);
             if (coeff.is_zero()) {
                 continue;
@@ -166,8 +267,7 @@ namespace opt {
 
 
     void model_based_opt::update_values(unsigned_vector const& bound_vars, unsigned_vector const& bound_trail) {
-        for (unsigned i = bound_trail.size(); i > 0; ) {
-            --i;
+        for (unsigned i = bound_trail.size(); i-- > 0; ) {
             unsigned x = bound_vars[i];
             row& r = m_rows[bound_trail[i]];
             rational val = r.m_coeff;
@@ -222,8 +322,7 @@ namespace opt {
         }        
         
         // update and check bounds for all other affected rows.
-        for (unsigned i = bound_trail.size(); i > 0; ) {
-            --i;
+        for (unsigned i = bound_trail.size(); i-- > 0; ) {
             unsigned x = bound_vars[i];
             unsigned_vector const& row_ids = m_var2row_ids[x];
             for (unsigned row_id : row_ids) {                
@@ -243,8 +342,7 @@ namespace opt {
         uint_set visited;
         m_above.reset();
         m_below.reset();
-        for (unsigned i = 0; i < row_ids.size(); ++i) {
-            unsigned row_id = row_ids[i];
+        for (unsigned row_id : row_ids) {
             SASSERT(row_id != m_objective_id);
             if (visited.contains(row_id)) {
                 continue;
@@ -291,8 +389,7 @@ namespace opt {
     rational model_based_opt::get_row_value(row const& r) const {
         vector<var> const& vars = r.m_vars;
         rational val = r.m_coeff;
-        for (unsigned i = 0; i < vars.size(); ++i) {
-            var const& v = vars[i];
+        for (var const& v : vars) {
             val += v.m_coeff * m_var2value[v.m_id];
         }
         return val;
@@ -300,14 +397,18 @@ namespace opt {
  
     rational model_based_opt::get_coefficient(unsigned row_id, unsigned var_id) const {
         row const& r = m_rows[row_id];
-        if (r.m_vars.empty()) {
+        return r.get_coefficient(var_id);
+    }
+
+    rational model_based_opt::row::get_coefficient(unsigned var_id) const {
+        if (m_vars.empty()) {
             return rational::zero();
         }
-        unsigned lo = 0, hi = r.m_vars.size();
+        unsigned lo = 0, hi = m_vars.size();
         while (lo < hi) {
             unsigned mid = lo + (hi - lo)/2;
             SASSERT(mid < hi);
-            unsigned id = r.m_vars[mid].m_id;
+            unsigned id = m_vars[mid].m_id;
             if (id == var_id) {
                 lo = mid;
                 break;
@@ -319,12 +420,12 @@ namespace opt {
                 hi = mid;
             }
         }
-        if (lo == r.m_vars.size()) {
+        if (lo == m_vars.size()) {
             return rational::zero();
         }
-        unsigned id = r.m_vars[lo].m_id;
+        unsigned id = m_vars[lo].m_id;
         if (id == var_id) {
-            return r.m_vars[lo].m_coeff;
+            return m_vars[lo].m_coeff;
         }
         else {
             return rational::zero();
@@ -386,7 +487,6 @@ namespace opt {
         SASSERT(a1 == get_coefficient(row_src, x));
         SASSERT(a1.is_pos());
         SASSERT(row_src != row_dst);                
-        SASSERT(m_rows[row_src].m_type == t_eq);
         if (!m_rows[row_dst].m_alive) return;
         rational a2 = get_coefficient(row_dst, x);
         mul(row_dst, a1);
@@ -593,40 +693,52 @@ namespace opt {
     }
     
     void model_based_opt::display(std::ostream& out) const {
-        for (unsigned i = 0; i < m_rows.size(); ++i) {
-            display(out, m_rows[i]);
+        for (auto const& r : m_rows) {
+            display(out, r);
         }
         for (unsigned i = 0; i < m_var2row_ids.size(); ++i) {
             unsigned_vector const& rows = m_var2row_ids[i];
             out << i << ": ";
-            for (unsigned j = 0; j < rows.size(); ++j) {
-                out << rows[j] << " ";
+            for (auto const& r : rows) {
+                out << r << " ";
             }
             out << "\n";
         }
     }        
 
-    void model_based_opt::display(std::ostream& out, row const& r) const {
-        vector<var> const& vars = r.m_vars;
-        out << (r.m_alive?"+":"-") << " ";
+    void model_based_opt::display(std::ostream& out, vector<var> const& vars, rational const& coeff) {
         for (unsigned i = 0; i < vars.size(); ++i) {
             if (i > 0 && vars[i].m_coeff.is_pos()) {
                 out << "+ ";
             }
             out << vars[i].m_coeff << "* v" << vars[i].m_id << " ";                
         }
-        if (r.m_coeff.is_pos()) {
-            out << " + " << r.m_coeff << " ";
+        if (coeff.is_pos()) {
+            out << " + " << coeff << " ";
         }
-        else if (r.m_coeff.is_neg()) {
-            out << r.m_coeff << " ";
-        }        
+        else if (coeff.is_neg()) {
+            out << coeff << " ";
+        }                
+    }
+
+    std::ostream& model_based_opt::display(std::ostream& out, row const& r) {
+        out << (r.m_alive?"+":"-") << " ";
+        display(out, r.m_vars, r.m_coeff);
         if (r.m_type == opt::t_mod) {
             out << r.m_type << " " << r.m_mod << " = 0; value: " << r.m_value  << "\n";
         }
         else {
             out << r.m_type << " 0; value: " << r.m_value  << "\n";
         }
+        return out;
+    }
+
+    std::ostream& model_based_opt::display(std::ostream& out, def const& r) {
+        display(out, r.m_vars, r.m_coeff);
+        if (!r.m_div.is_one()) {
+            out << " / " << r.m_div;
+        }
+        return out;
     }
 
     unsigned model_based_opt::add_var(rational const& value, bool is_int) {
@@ -648,10 +760,10 @@ namespace opt {
         r.m_vars.append(coeffs.size(), coeffs.c_ptr());
         bool is_int_row = true;
         std::sort(r.m_vars.begin(), r.m_vars.end(), var::compare());
-        for (unsigned i = 0; i < coeffs.size(); ++i) {
-            val += m_var2value[coeffs[i].m_id] * coeffs[i].m_coeff;
-            SASSERT(!is_int(coeffs[i].m_id) || coeffs[i].m_coeff.is_int());
-            is_int_row &= is_int(coeffs[i].m_id);
+        for (auto const& c : coeffs) {
+            val += m_var2value[c.m_id] * c.m_coeff;
+            SASSERT(!is_int(c.m_id) || c.m_coeff.is_int());
+            is_int_row &= is_int(c.m_id);
         }
         r.m_alive = true;
         r.m_coeff = c;
@@ -738,7 +850,7 @@ namespace opt {
     //       t0 <= s for each s (M inequalities).
     // If N >= M the construction is symmetric.
     // 
-    void model_based_opt::project(unsigned x) {
+    model_based_opt::def model_based_opt::project(unsigned x, bool compute_def) {
         unsigned_vector& lub_rows = m_lub;
         unsigned_vector& glb_rows = m_glb;
         unsigned_vector& mod_rows = m_mod;
@@ -802,28 +914,41 @@ namespace opt {
         }
 
         if (!mod_rows.empty()) {
-            solve_mod(x, mod_rows);
-            return;
+            return solve_mod(x, mod_rows, compute_def);
         }
 
         if (eq_row != UINT_MAX) {
-            solve_for(eq_row, x);
-            return;
+            return solve_for(eq_row, x, compute_def);
         }
-        
+
+        def result;        
         unsigned lub_size = lub_rows.size();
         unsigned glb_size = glb_rows.size();
         unsigned row_index = (lub_size <= glb_size) ? lub_index : glb_index;
-        glb_rows.append(lub_rows);  
         
         // There are only upper or only lower bounds.
         if (row_index == UINT_MAX) {
-            for (unsigned row_id : glb_rows) {
-                SASSERT(m_rows[row_id].m_alive);
-                SASSERT(!get_coefficient(row_id, x).is_zero());
-                retire_row(row_id);
+            if (compute_def) {
+                if (lub_index != UINT_MAX) {
+                    result = solve_for(lub_index, x, true);
+                }
+                else if (glb_index != UINT_MAX) {
+                    result = solve_for(glb_index, x, true);                
+                }
             }
-            return;
+            else {
+                for (unsigned row_id : lub_rows) retire_row(row_id);
+                for (unsigned row_id : glb_rows) retire_row(row_id);
+            }
+            return result;
+        }
+
+        SASSERT(lub_index != UINT_MAX);
+        SASSERT(glb_index != UINT_MAX);
+        if (compute_def) {
+            def d1(m_rows[lub_index], x);
+            def d2(m_rows[lub_index], x);
+            result = (d1 + d2)/2;
         }
 
         // The number of matching lower and upper bounds is small.
@@ -832,10 +957,9 @@ namespace opt {
             (!is_int(x) || lub_is_unit || glb_is_unit)) {
             for (unsigned i = 0; i < lub_size; ++i) {
                 unsigned row_id1 = lub_rows[i];
-                bool last = i + 1 == lub_rows.size();
+                bool last = i + 1 == lub_size;
                 rational coeff = get_coefficient(row_id1, x);
-                for (unsigned j = 0; j < glb_size; ++j) {
-                    unsigned row_id2 = glb_rows[j];
+                for (unsigned row_id2 : glb_rows) {
                     if (last) {
                         resolve(row_id1, coeff, row_id2, x);
                     }
@@ -845,20 +969,25 @@ namespace opt {
                     }
                 }
             }
-            for (unsigned i = 0; i < lub_size; ++i) {
-                retire_row(lub_rows[i]);
-            }
-            return;
+            for (unsigned row_id : lub_rows) retire_row(row_id);
+
+            return result;
         }
 
         // General case.
         rational coeff = get_coefficient(row_index, x);
+        for (unsigned row_id : lub_rows) {
+            if (row_id != row_index) {
+                resolve(row_index, coeff, row_id, x);
+            }
+        }
         for (unsigned row_id : glb_rows) {
             if (row_id != row_index) {
                 resolve(row_index, coeff, row_id, x);
             }
         }
         retire_row(row_index);
+        return result;
     }
 
     // 
@@ -876,7 +1005,7 @@ namespace opt {
     // x := D*x' + u
     // 
 
-    void model_based_opt::solve_mod(unsigned x, unsigned_vector const& mod_rows) {
+    model_based_opt::def model_based_opt::solve_mod(unsigned x, unsigned_vector const& mod_rows, bool compute_def) {
         SASSERT(!mod_rows.empty());
         rational D(1);
         for (unsigned idx : mod_rows) {
@@ -914,7 +1043,11 @@ namespace opt {
                 visited.insert(row_id);
             }
         }
-        project(y);
+        def result = project(y, compute_def);
+        if (compute_def) {
+            result = (result * D) + u;
+        }
+        return result;
     }
 
     // update row with: x |-> C
@@ -961,16 +1094,21 @@ namespace opt {
     // 3x + t = 0 & 7 | (c*x + s) & ax <= u 
     // 3 | -t  & 21 | (-ct + 3s) & a-t <= 3u
 
-    void model_based_opt::solve_for(unsigned row_id1, unsigned x) {
+    model_based_opt::def model_based_opt::solve_for(unsigned row_id1, unsigned x, bool compute_def) {
+        TRACE("opt", tout << "v" << x << "\n" << m_rows[row_id1] << "\n";);
         rational a = get_coefficient(row_id1, x), b;
+        ineq_type ty = m_rows[row_id1].m_type;
         SASSERT(!a.is_zero());
-        SASSERT(m_rows[row_id1].m_type == t_eq);
         SASSERT(m_rows[row_id1].m_alive);
         if (a.is_neg()) {
             a.neg();
             m_rows[row_id1].neg();
         }
         SASSERT(a.is_pos());
+        if (ty == t_lt) {
+            SASSERT(compute_def);
+            m_rows[row_id1].m_coeff += a;            
+        }
         if (m_var2is_int[x] && !a.is_one()) {
             row& r1 = m_rows[row_id1];
             vector<var> coeffs;
@@ -983,8 +1121,7 @@ namespace opt {
         visited.insert(row_id1);
         for (unsigned row_id2 : row_ids) {
             if (!visited.contains(row_id2)) {                
-                visited.insert(row_id2);
-                
+                visited.insert(row_id2);                
                 b = get_coefficient(row_id2, x);
                 if (!b.is_zero()) {
                     row& dst = m_rows[row_id2];
@@ -1002,14 +1139,19 @@ namespace opt {
                 }
             }
         }
+        // TBD: -t div a
+        def result(m_rows[row_id1], x);
         retire_row(row_id1);
+        return result;
     }
-
-    void model_based_opt::project(unsigned num_vars, unsigned const* vars) {
+    
+    vector<model_based_opt::def> model_based_opt::project(unsigned num_vars, unsigned const* vars, bool compute_def) {
+        vector<def> result;
         for (unsigned i = 0; i < num_vars; ++i) {
-            project(vars[i]);
+            result.push_back(project(vars[i], compute_def));
             TRACE("opt", display(tout << "After projecting: v" << vars[i] << "\n"););
         }
+        return result;
     }
 
 }
diff --git a/src/math/simplex/model_based_opt.h b/src/math/simplex/model_based_opt.h
index 9546529f2..268d3d81d 100644
--- a/src/math/simplex/model_based_opt.h
+++ b/src/math/simplex/model_based_opt.h
@@ -60,6 +60,23 @@ namespace opt {
             void reset() { m_vars.reset(); m_coeff.reset(); m_value.reset(); }
 
             void neg() { for (var & v : m_vars) v.m_coeff.neg(); m_coeff.neg(); m_value.neg(); }
+            rational get_coefficient(unsigned x) const;
+        };
+
+        // A definition is a linear term of the form (vars + coeff) / div 
+        struct def {
+            def(): m_div(1) {}
+            def(row const& r, unsigned x);
+            def(def const& other): m_vars(other.m_vars), m_coeff(other.m_coeff), m_div(other.m_div) {}
+            vector<var> m_vars;
+            rational    m_coeff;
+            rational    m_div; 
+            def operator+(def const& other) const;
+            def operator/(unsigned n) const { return *this / rational(n); }
+            def operator/(rational const& n) const;
+            def operator*(rational const& n) const;
+            def operator+(rational const& n) const;
+            void normalize();
         };
 
     private:
@@ -121,12 +138,12 @@ namespace opt {
 
         void update_value(unsigned x, rational const& val);
 
-        void project(unsigned var);
+        def project(unsigned var, bool compute_def);
 
-        void solve_for(unsigned row_id, unsigned x);
-
-        void solve_mod(unsigned x, unsigned_vector const& mod_rows);
+        def solve_for(unsigned row_id, unsigned x, bool compute_def);
 
+        def solve_mod(unsigned x, unsigned_vector const& mod_rows, bool compute_def);
+        
         bool is_int(unsigned x) const { return m_var2is_int[x]; }
 
         void retire_row(unsigned row_id);
@@ -163,7 +180,7 @@ namespace opt {
         // 
         // Project set of variables from inequalities.
         //
-        void project(unsigned num_vars, unsigned const* vars);
+        vector<def> project(unsigned num_vars, unsigned const* vars, bool compute_def);
 
         //
         // Extract current rows (after projection).
@@ -171,13 +188,17 @@ namespace opt {
         void get_live_rows(vector<row>& rows);
 
         void display(std::ostream& out) const;
-        void display(std::ostream& out, row const& r) const;
+        static std::ostream& display(std::ostream& out, row const& r);
+        static std::ostream& display(std::ostream& out, def const& r);
+        static void display(std::ostream& out, vector<var> const& vars, rational const& coeff);
 
     };
 
 }
 
 std::ostream& operator<<(std::ostream& out, opt::ineq_type ie);
+inline std::ostream& operator<<(std::ostream& out, opt::model_based_opt::def const& d) { return opt::model_based_opt::display(out, d); }
+inline std::ostream& operator<<(std::ostream& out, opt::model_based_opt::row const& r) { return opt::model_based_opt::display(out, r); }
 
 inline std::ostream& operator<<(std::ostream& out, opt::model_based_opt::var const v) { return out << "v" << v.m_id; }
 
diff --git a/src/qe/qe_arith.cpp b/src/qe/qe_arith.cpp
index 0623bc8bb..0b5ad72ca 100644
--- a/src/qe/qe_arith.cpp
+++ b/src/qe/qe_arith.cpp
@@ -283,14 +283,29 @@ namespace qe {
         typedef opt::model_based_opt::row row;
         typedef vector<var> vars;
 
+        vector<def> project(model& model, app_ref_vector& vars, expr_ref_vector& lits) {
+            return project(model, vars, lits, true);
+        }
+
         void operator()(model& model, app_ref_vector& vars, expr_ref_vector& fmls) {
+            project(model, vars, fmls, false);
+        }
+
+        expr_ref var2expr(ptr_vector<expr> const& index2expr, var const& v) {
+            expr_ref t(index2expr[v.m_id], m);
+            if (!v.m_coeff.is_one()) {
+                t = a.mk_mul(a.mk_numeral(v.m_coeff, a.is_int(t)), t);
+            }
+            return t;
+        }
+
+        vector<def> project(model& model, app_ref_vector& vars, expr_ref_vector& fmls, bool compute_def) {
             bool has_arith = false;
-            for (unsigned i = 0; !has_arith && i < vars.size(); ++i) {
-                expr* v = vars[i].get();
+            for (expr* v : vars) {
                 has_arith |= is_arith(v);
             }
             if (!has_arith) {
-                return;
+                return vector<def>();
             }
             model_evaluator eval(model);
             // eval.set_model_completion(true);
@@ -359,7 +374,7 @@ namespace qe {
                       tout << "v" << v << " " << mk_pp(index2expr[v], m) << "\n";
                   }
                   mbo.display(tout););
-            mbo.project(real_vars.size(), real_vars.c_ptr());
+            vector<opt::model_based_opt::def> defs = mbo.project(real_vars.size(), real_vars.c_ptr(), compute_def);
             TRACE("qe", mbo.display(tout););
             vector<row> rows;
             mbo.get_live_rows(rows);
@@ -419,6 +434,30 @@ namespace qe {
                 val = eval(t);
                 CTRACE("qe", !m.is_true(val), tout << "Evaluated " << t << " to " << val << "\n";);
             }
+            vector<def> result;
+            if (compute_def) {
+                SASSERT(defs.size() == real_vars.size());
+                for (unsigned i = 0; i < defs.size(); ++i) {
+                    auto const& d = defs[i];
+                    expr* x = index2expr[real_vars[i]];
+                    bool is_int = a.is_int(x);
+                    expr_ref_vector ts(m);
+                    expr_ref t(m);
+                    for (var const& v : d.m_vars) {
+                        ts.push_back(var2expr(index2expr, v));
+                    }
+                    ts.push_back(a.mk_numeral(d.m_coeff, is_int));
+                    t = a.mk_add(ts.size(), ts.c_ptr());
+                    if (!d.m_div.is_one() && is_int) {
+                        t = a.mk_idiv(t, a.mk_numeral(d.m_div, is_int));
+                    }
+                    else if (!d.m_div.is_one() && !is_int) {
+                        t = a.mk_div(t, a.mk_numeral(d.m_div, is_int));
+                    }
+                    result.push_back(def(expr_ref(x, m), t));
+                }
+            }
+            return result;
         }        
 
         opt::inf_eps maximize(expr_ref_vector const& fmls0, model& mdl, app* t, expr_ref& ge, expr_ref& gt) {
@@ -534,6 +573,10 @@ namespace qe {
         (*m_imp)(model, vars, lits);
     }
 
+    vector<def> arith_project_plugin::project(model& model, app_ref_vector& vars, expr_ref_vector& lits) {
+        return m_imp->project(model, vars, lits);
+    }
+
     bool arith_project_plugin::solve(model& model, app_ref_vector& vars, expr_ref_vector& lits) {
         return m_imp->solve(model, vars, lits);
     }
diff --git a/src/qe/qe_arith.h b/src/qe/qe_arith.h
index 8e2400ac7..c01d7bbb6 100644
--- a/src/qe/qe_arith.h
+++ b/src/qe/qe_arith.h
@@ -30,7 +30,7 @@ namespace qe {
         bool solve(model& model, app_ref_vector& vars, expr_ref_vector& lits) override;
         family_id get_family_id() override;
         void operator()(model& model, app_ref_vector& vars, expr_ref_vector& lits) override;
-
+        vector<def> project(model& model, app_ref_vector& vars, expr_ref_vector& lits) override;
 
         opt::inf_eps maximize(expr_ref_vector const& fmls, model& mdl, app* t, expr_ref& ge, expr_ref& gt);
     };
diff --git a/src/qe/qe_arrays.cpp b/src/qe/qe_arrays.cpp
index 651c167b6..138aed1df 100644
--- a/src/qe/qe_arrays.cpp
+++ b/src/qe/qe_arrays.cpp
@@ -1419,7 +1419,9 @@ namespace qe {
               );
     }
 
-
+    vector<def> array_project_plugin::project(model& model, app_ref_vector& vars, expr_ref_vector& lits) {
+        return vector<def>();
+    }
 
 };
 
diff --git a/src/qe/qe_arrays.h b/src/qe/qe_arrays.h
index f2d69b564..3bb90335d 100644
--- a/src/qe/qe_arrays.h
+++ b/src/qe/qe_arrays.h
@@ -36,6 +36,7 @@ namespace qe {
         bool solve(model& model, app_ref_vector& vars, expr_ref_vector& lits) override;
         void operator()(model& model, app_ref_vector& vars, expr_ref& fml, app_ref_vector& aux_vars, bool reduce_all_selects);
         family_id get_family_id() override;
+        vector<def> project(model& model, app_ref_vector& vars, expr_ref_vector& lits) override;
     };
 
 };
diff --git a/src/qe/qe_datatypes.cpp b/src/qe/qe_datatypes.cpp
index 770beb59d..87ae97857 100644
--- a/src/qe/qe_datatypes.cpp
+++ b/src/qe/qe_datatypes.cpp
@@ -300,6 +300,9 @@ namespace qe {
         return m_imp->solve(model, vars, lits);
     }
 
+    vector<def> datatype_project_plugin::project(model& model, app_ref_vector& vars, expr_ref_vector& lits) {
+        return vector<def>();
+    }
     
     family_id datatype_project_plugin::get_family_id() {
         return m_imp->dt.get_family_id();
diff --git a/src/qe/qe_datatypes.h b/src/qe/qe_datatypes.h
index 9c8eb5a76..0483f4cce 100644
--- a/src/qe/qe_datatypes.h
+++ b/src/qe/qe_datatypes.h
@@ -35,6 +35,7 @@ namespace qe {
         bool operator()(model& model, app* var, app_ref_vector& vars, expr_ref_vector& lits) override;
         bool solve(model& model, app_ref_vector& vars, expr_ref_vector& lits) override;
         family_id get_family_id() override;
+        vector<def> project(model& model, app_ref_vector& vars, expr_ref_vector& lits) override;
     };
 
 };
diff --git a/src/qe/qe_mbp.h b/src/qe/qe_mbp.h
index bddec8065..c35351fb3 100644
--- a/src/qe/qe_mbp.h
+++ b/src/qe/qe_mbp.h
@@ -31,14 +31,32 @@ namespace qe {
 
     struct cant_project {};
 
+    struct def {
+        expr_ref var, term;
+        def(expr_ref& v, expr_ref& t): var(v), term(t) {}
+    };
+
     class project_plugin {
     public:
         virtual ~project_plugin() {}
         virtual bool operator()(model& model, app* var, app_ref_vector& vars, expr_ref_vector& lits) = 0;
+        /**
+           \brief partial solver.
+        */
         virtual bool solve(model& model, app_ref_vector& vars, expr_ref_vector& lits) = 0;
         virtual family_id get_family_id() = 0;
+
         virtual void operator()(model& model, app_ref_vector& vars, expr_ref_vector& lits) { };
 
+        /**
+           \brief project vars modulo model, return set of definitions for eliminated variables.
+           - vars in/out: returns variables that were not eliminated
+           - lits in/out: returns projected literals
+           - returns set of definitions
+             (TBD: in triangular form, the last definition can be substituted into definitions that come before)
+        */
+        virtual vector<def> project(model& model, app_ref_vector& vars, expr_ref_vector& lits) = 0;
+
         static expr_ref pick_equality(ast_manager& m, model& model, expr* t);
         static void erase(expr_ref_vector& lits, unsigned& i);
         static void push_back(expr_ref_vector& lits, expr* lit);
diff --git a/src/test/CMakeLists.txt b/src/test/CMakeLists.txt
index 1dc0290d8..7367a52ff 100644
--- a/src/test/CMakeLists.txt
+++ b/src/test/CMakeLists.txt
@@ -81,7 +81,6 @@ add_executable(test-z3
   optional.cpp
   parray.cpp
   pb2bv.cpp
-  pdr.cpp
   permutation.cpp
   polynomial.cpp
   polynorm.cpp
diff --git a/src/test/arith_rewriter.cpp b/src/test/arith_rewriter.cpp
index a279d58b2..d0a110c4f 100644
--- a/src/test/arith_rewriter.cpp
+++ b/src/test/arith_rewriter.cpp
@@ -10,7 +10,6 @@ Copyright (c) 2015 Microsoft Corporation
 #include "ast/reg_decl_plugins.h"
 #include "ast/rewriter/th_rewriter.h"
 #include "model/model.h"
-#include "muz/pdr/pdr_util.h"
 #include "parsers/smt2/smt2parser.h"
 
 
@@ -53,13 +52,9 @@ void tst_arith_rewriter() {
     expr_ref fml = parse_fml(m, example1);
     rw(fml);
     std::cout << mk_pp(fml, m) << "\n";
-    pdr::normalize_arithmetic(fml);
-    std::cout << mk_pp(fml, m) << "\n";
 
 
     fml = parse_fml(m, example2);
     rw(fml);
     std::cout << mk_pp(fml, m) << "\n";
-    pdr::normalize_arithmetic(fml);
-    std::cout << mk_pp(fml, m) << "\n";
 }
diff --git a/src/test/main.cpp b/src/test/main.cpp
index 41f051f24..d330610d9 100644
--- a/src/test/main.cpp
+++ b/src/test/main.cpp
@@ -239,7 +239,6 @@ int main(int argc, char ** argv) {
     TST(theory_pb);
     TST(simplex);
     TST(sat_user_scope);
-    TST(pdr);
     TST_ARGV(ddnf);
     TST(ddnf1);
     TST(model_evaluator);
diff --git a/src/test/model_based_opt.cpp b/src/test/model_based_opt.cpp
index 7b3b95afd..f2f6004ff 100644
--- a/src/test/model_based_opt.cpp
+++ b/src/test/model_based_opt.cpp
@@ -218,6 +218,12 @@ static void test4() {
     std::cout << "u: " << mbo.get_value(u) << "\n";
 }
 
+static void display_project(vector<opt::model_based_opt::def> const& defs) {
+    for (auto const& d : defs) {
+        std::cout << d << "\n";
+    }
+}
+
 static void test5() {
     opt::model_based_opt mbo;
     unsigned x = mbo.add_var(rational(2));
@@ -231,13 +237,13 @@ static void test5() {
     add_ineq(mbo, z, 1, u, -1, 1, opt::t_le);
 
     unsigned vars[2] = { y, z };
-    mbo.project(1, vars);    
+    display_project(mbo.project(1, vars, true));    
     mbo.display(std::cout);
 
-    mbo.project(1, vars);    
+    display_project(mbo.project(1, vars, true));    
     mbo.display(std::cout);
 
-    mbo.project(1, vars+1);
+    display_project(mbo.project(1, vars+1, true));
     mbo.display(std::cout);
 
     vector<opt::model_based_opt::row> rows;
@@ -263,7 +269,7 @@ static void test6() {
     add_ineq(mbo, y, 1, w, -1, 1, opt::t_le);
 
     mbo.display(std::cout);
-    mbo.project(1, &y);
+    display_project(mbo.project(1, &y, true));
     mbo.display(std::cout);
 }
 
@@ -285,7 +291,7 @@ static void test7() {
     add_ineq(mbo, y, 1, w, -1, 1, opt::t_lt);
 
     mbo.display(std::cout);
-    mbo.project(1, &y);
+    display_project(mbo.project(1, &y, true));
     mbo.display(std::cout);
 }
 
@@ -306,7 +312,7 @@ static void test8() {
     add_ineq(mbo, y, 1, v, -1, 1, opt::t_le);
 
     mbo.display(std::cout);
-    mbo.project(1, &y);
+    display_project(mbo.project(1, &y, true));
     mbo.display(std::cout);
 }
 
@@ -327,7 +333,7 @@ static void test9() {
     add_ineq(mbo, y, 1, v, -1, 1, opt::t_le);
 
     mbo.display(std::cout);
-    mbo.project(1, &y);
+    display_project(mbo.project(1, &y, true));
     mbo.display(std::cout);
 }
 
@@ -348,9 +354,9 @@ static void test10() {
     add_ineq(mbo, y, 3, v, -6, 1, opt::t_le);
 
     mbo.display(std::cout);
-    mbo.project(1, &y);
+    display_project(mbo.project(1, &y, true));
     mbo.display(std::cout);
-    mbo.project(1, &x0);
+    display_project(mbo.project(1, &x0, true));
     mbo.display(std::cout);
 }
 
@@ -358,7 +364,6 @@ static void test10() {
 
 void tst_model_based_opt() {
     test10();
-    return;
     check_random_ineqs();
     test1();
     test2();
diff --git a/src/test/pdr.cpp b/src/test/pdr.cpp
deleted file mode 100644
index 45d9eea7e..000000000
--- a/src/test/pdr.cpp
+++ /dev/null
@@ -1,128 +0,0 @@
-
-/*++
-Copyright (c) 2015 Microsoft Corporation
-
---*/
-
-#include "muz/pdr/pdr_context.h"
-#include "ast/reg_decl_plugins.h"
-
-
-using namespace pdr;
-
-static expr_ref mk_state(expr_ref_vector const& states, random_gen& rand) {
-    expr_ref result(states.get_manager());
-    result = states[rand(states.size())];
-    return result;
-}
-
-
-struct test_model_search {
-    struct init_test {
-        init_test(func_decl_ref& fn) {
-            ast_manager& m = fn.get_manager();
-            reg_decl_plugins(m);
-            fn = m.mk_const_decl(symbol("f"), m.mk_bool_sort());
-        }
-    };
-
-    ast_manager       m;
-    smt_params        m_smt_params;
-    fixedpoint_params fp_params;
-    context           ctx;
-    manager           pm;
-    func_decl_ref     fn;
-    init_test         initt;
-    pred_transformer  pt;
-    random_gen        rand;
-    model_search      search;
-    expr_ref_vector   states;  
-
-
-    test_model_search():
-        ctx(m_smt_params, fp_params, m),
-        pm(m_smt_params, 10, m),
-        fn(m),
-        initt(fn),
-        pt(ctx, pm, fn),
-        rand(10),
-        search(true),
-        states(m) {
-    }
-
-    void add_tree(model_node* parent, bool force_goal) {
-        unsigned level = parent->level();
-        search.add_leaf(*parent);
-        expr_ref state(m);
-        if (level > 0 && (force_goal || parent->is_goal())) {
-            search.remove_goal(*parent);
-
-            state = mk_state(states, rand);
-            add_tree(alloc(model_node, parent, state, pt, level-1), false);
-
-            state = mk_state(states, rand);
-            add_tree(alloc(model_node, parent, state, pt, level-1), false);
-            parent->check_pre_closed();
-        }
-    }
-
-    bool mutate() {
-        model_node* leaf = search.next();
-        if (!leaf) return false;
-        unsigned level = leaf->level();
-        if (level == 0) {
-            if (rand(2) == 0) {
-                leaf->display(std::cout << "backtrack to grandparent\n", 1);
-                search.backtrack_level(false, *leaf->parent());
-            }
-            else {
-                leaf->display(std::cout << "backtrack to parent\n", 1);
-                search.backtrack_level(false, *leaf);
-            }
-        }
-        else {
-            leaf->display(std::cout << "grow tree\n", 1);
-            add_tree(leaf, true);
-        }
-        return true;
-    }
-
-    void init() {
-        std::cout << "pdr state-hash search tree\n";
-        
-        expr_ref state(m);
-        func_decl_ref fn(m);
-        for (unsigned i = 0; i < 10; ++i) {
-            std::ostringstream strm;
-            strm << "s" << i;
-            state = m.mk_const(symbol(strm.str().c_str()), m.mk_bool_sort());
-            fn = to_app(state)->get_decl();
-            states.push_back(state);
-        }
-        state = states[0].get();
-        unsigned level = 4;
-        for(unsigned n = 0; n < 100; ++n) {
-            state = mk_state(states, rand);
-            model_node* root = alloc(model_node, nullptr, state, pt, level);
-            search.set_root(root);
-            add_tree(root, false);
-            search.display(std::cout);
-            
-            while (true) {
-                search.well_formed();
-                if (!mutate()) break;
-                search.display(std::cout);
-            }
-            search.reset();
-            //++level;
-        }
-        search.reset();        
-    }
-
-};
-
-void tst_pdr() {
-    test_model_search test;
-
-    test.init();
-}