z3-src 415.4.3

/*++
Copyright (c) 2006 Microsoft Corporation

Module Name:

    theory_dense_diff_logic_def.h

Abstract:

    <abstract>

Author:

    Leonardo de Moura (leonardo) 2008-05-16.

Revision History:

--*/
#pragma once

#include "smt/smt_context.h"
#include "smt/theory_dense_diff_logic.h"
#include "ast/ast_pp.h"
#include "smt/smt_model_generator.h"
#include "math/simplex/simplex.h"
#include "math/simplex/simplex_def.h"

namespace smt {

    template<typename Ext>
    theory_dense_diff_logic<Ext>::theory_dense_diff_logic(context& ctx):
        theory(ctx, ctx.get_manager().mk_family_id("arith")),
        m_params(ctx.get_fparams()),
        m_autil(ctx.get_manager()),
        m_arith_eq_adapter(*this, m_autil),
        m_non_diff_logic_exprs(false),
        m_var_value_table(DEFAULT_HASHTABLE_INITIAL_CAPACITY, var_value_hash(*this), var_value_eq(*this)) {
        m_edges.push_back(edge());
    }

    template<typename Ext>
    theory* theory_dense_diff_logic<Ext>::mk_fresh(context * new_ctx) { 
        return alloc(theory_dense_diff_logic<Ext>, *new_ctx);
    }

    template<typename Ext>
    inline app * theory_dense_diff_logic<Ext>::mk_zero_for(expr * n) {
        return m_autil.mk_numeral(rational(0), n->get_sort());
    }

    template<typename Ext>
    theory_var theory_dense_diff_logic<Ext>::mk_var(enode * n) {
        theory_var v = theory::mk_var(n);
        bool is_int  = m_autil.is_int(n->get_expr());
        m_is_int.push_back(is_int);
        m_f_targets.push_back(f_target());
        for (auto& rows : m_matrix) {
            rows.push_back(cell());
        }
        m_matrix.push_back(row());
        row & r = m_matrix.back();
        SASSERT(r.empty());
        r.resize(v+1);
        cell & c    = m_matrix[v][v];
        c.m_edge_id = self_edge_id;
        c.m_distance.reset();
        SASSERT(check_vector_sizes());
        ctx.attach_th_var(n, this, v);
        return v;
    }

    template<typename Ext>
    theory_var theory_dense_diff_logic<Ext>::internalize_term_core(app * n) {
        if (ctx.e_internalized(n)) {
            enode * e    = ctx.get_enode(n);
            if (is_attached_to_var(e))
                return e->get_th_var(get_id());
        }

        rational _k;
        if (m_autil.is_add(n) && to_app(n)->get_num_args() == 2 && m_autil.is_numeral(to_app(n)->get_arg(0), _k)) {
            numeral k(_k);
            if (m_params.m_arith_reflect)
                internalize_term_core(to_app(to_app(n)->get_arg(0)));
            theory_var s = internalize_term_core(to_app(to_app(n)->get_arg(1)));
            if (null_theory_var == s)
                return s;
            enode * e    = ctx.mk_enode(n, !m_params.m_arith_reflect, false, true);
            theory_var v = mk_var(e);
            add_edge(s, v, k, null_literal);
            k.neg();
            add_edge(v, s, k, null_literal);
            return v;
        }
        else if (m_autil.is_numeral(n, _k)) {
            enode * e    = ctx.mk_enode(n, false, false, true);
            // internalizer is marking enodes as interpreted whenever the associated ast is a value and a constant.
            // e->mark_as_interpreted();
            theory_var v = mk_var(e);
            if (_k.is_zero()) 
                return v;
            theory_var z = internalize_term_core(mk_zero_for(n));
            numeral k(_k);
            add_edge(z, v, k, null_literal);
            k.neg();
            add_edge(v, z, k, null_literal);
            return v;
        }
        else if (!m_autil.is_arith_expr(n)) {
            if (!ctx.e_internalized(n))
                ctx.internalize(n, false);
            enode * e    = ctx.get_enode(n);
            if (!is_attached_to_var(e))
                return mk_var(e);
            else
                return e->get_th_var(get_id());
        }
        else {
            return null_theory_var;
        }
    }

    template<typename Ext>
    void theory_dense_diff_logic<Ext>::found_non_diff_logic_expr(expr * n) {
        if (!m_non_diff_logic_exprs) {
            TRACE(non_diff_logic, tout << "found non diff logic expression:\n" << mk_pp(n, m) << "\n";);
            ctx.push_trail(value_trail<bool>(m_non_diff_logic_exprs));
            IF_VERBOSE(0, verbose_stream() << "(smt.diff_logic: non-diff logic expression " << mk_pp(n, m) << ")\n";); 
            m_non_diff_logic_exprs = true;
        }
    }
                          
    template<typename Ext>
    bool theory_dense_diff_logic<Ext>::internalize_atom(app * n, bool gate_ctx) {
        if (memory::above_high_watermark()) {
            found_non_diff_logic_expr(n); // little hack... TODO: change to no_memory and return l_undef if SAT
            return false;
        }
        TRACE(ddl, tout << "internalizing atom:\n" << mk_pp(n, m) << "\n";);
        SASSERT(!ctx.b_internalized(n));
        SASSERT(m_autil.is_le(n) || m_autil.is_ge(n));
        theory_var source, target;
        SASSERT(m_autil.is_le(n) || m_autil.is_ge(n));
        app * lhs      = to_app(n->get_arg(0));
        app * rhs      = to_app(n->get_arg(1));
        if (!m_autil.is_numeral(rhs)) {
            found_non_diff_logic_expr(n);
            return false;
        }
        rational _k;
        m_autil.is_numeral(rhs, _k);
        numeral offset(_k);
        app * s, * t;
        expr *arg1, *arg2;
        if (m_autil.is_add(lhs, arg1, arg2) && is_times_minus_one(arg2, s) && !m_autil.is_arith_expr(s) && !m_autil.is_arith_expr(arg1)) {
            t = to_app(arg1);
        }
        else if (m_autil.is_add(lhs, arg1, arg2) && is_times_minus_one(arg1, s) && !m_autil.is_arith_expr(s) && !m_autil.is_arith_expr(arg2)) {
            t = to_app(arg2);
        }
        else if (m_autil.is_mul(lhs, arg1, arg2) && m_autil.is_minus_one(arg1) && !m_autil.is_arith_expr(arg2)) {
            s = to_app(arg2);
            t = mk_zero_for(s);
        }
        else if (!m_autil.is_arith_expr(lhs)) {
            t = to_app(lhs);
            s = mk_zero_for(t);
        }
        else {
            TRACE(ddl, tout << "failed to internalize:\n" << mk_pp(n, m) << "\n";);
            found_non_diff_logic_expr(n);
            return false;
        }
        TRACE(arith, tout << expr_ref(lhs, m) << " " << expr_ref(s, m) << " " << expr_ref(t, m) << "\n";);
        source = internalize_term_core(s);
        target = internalize_term_core(t);
        if (source == null_theory_var || target == null_theory_var) {
            TRACE(ddl, tout << "failed to internalize:\n" << mk_pp(n, m) << "\n";);
            found_non_diff_logic_expr(n);
            return false;
        }
        SASSERT(source != null_theory_var && target != null_theory_var);
        if (m_autil.is_ge(n)) {
            std::swap(source, target);
            offset.neg();
        }
		if (ctx.b_internalized(n)) return true;
        bool_var bv = ctx.mk_bool_var(n);
        ctx.set_var_theory(bv, get_id());
        atom * a    = alloc(atom, bv, source, target, offset);
        m_atoms.push_back(a);
        m_bv2atoms.setx(bv, a, 0);
        m_matrix[source][target].m_occs.push_back(a);
        m_matrix[target][source].m_occs.push_back(a);
        TRACE(ddl, tout << "succeeded internalizing:\n" << mk_pp(n, m) << "\n";);
        return true;
    }

    template<typename Ext>
    void theory_dense_diff_logic<Ext>::mk_clause(literal l1, literal l2) {
        ctx.mk_th_axiom(get_id(), l1, l2);
    }

    template<typename Ext>
    void theory_dense_diff_logic<Ext>::mk_clause(literal l1, literal l2, literal l3) {
        ctx.mk_th_axiom(get_id(), l1, l2, l3);
    }

    template<typename Ext>
    bool theory_dense_diff_logic<Ext>::internalize_term(app * term) {        
        if (memory::above_high_watermark()) {
            found_non_diff_logic_expr(term); // little hack... TODO: change to no_memory and return l_undef if SAT
            return false;
        }
        TRACE(ddl, tout << "internalizing term: " << mk_pp(term, m) << "\n";);
        theory_var v = internalize_term_core(term);
        TRACE(ddl, tout << mk_pp(term, m) << "\ninternalization result: " << (v != null_theory_var) << "\n";);
        if (v == null_theory_var)
            found_non_diff_logic_expr(term);
        return v != null_theory_var;
    }

    template<typename Ext>
    void theory_dense_diff_logic<Ext>::internalize_eq_eh(app * atom, bool_var v) {
        TRACE(ddl, tout << "eq-eh: " << mk_pp(atom, m) << "\n";);
        if (memory::above_high_watermark())
            return;
        app * lhs      = to_app(atom->get_arg(0));
        app * rhs      = to_app(atom->get_arg(1));
        app * s;
        if (m_autil.is_add(lhs) && to_app(lhs)->get_num_args() == 2 && is_times_minus_one(to_app(lhs)->get_arg(1), s) 
            && m_autil.is_numeral(rhs)) {
            // force axioms for (= (+ x (* -1 y)) k)
            // this is necessary because (+ x (* -1 y)) is not a diff logic term.
            m_arith_eq_adapter.mk_axioms(ctx.get_enode(lhs), ctx.get_enode(rhs));
            return;
        }

        if (m_params.m_arith_eager_eq_axioms) {
            enode * n1 = ctx.get_enode(lhs);
            enode * n2 = ctx.get_enode(rhs);
            if (n1->get_th_var(get_id()) != null_theory_var &&
                n2->get_th_var(get_id()) != null_theory_var) {
                m_arith_eq_adapter.mk_axioms(n1, n2); 
                return;
            }
        }
    }
    
    template<typename Ext>
    void theory_dense_diff_logic<Ext>::apply_sort_cnstr(enode * n, sort * s) {
        // do nothing...
    }

    template<typename Ext>
    void theory_dense_diff_logic<Ext>::assign_eh(bool_var v, bool is_true) {
        if (ctx.has_th_justification(v, get_id())) {
            TRACE(ddl, tout << "ignoring atom propagated by the theory.\n";);
            return;
        }
        atom * a = m_bv2atoms.get(v, 0);
        if (!a) {
            SASSERT(m.is_eq(ctx.bool_var2expr(v)));
            return;
        }
        m_stats.m_num_assertions++;
        literal l    = literal(v, !is_true);
        theory_var s = a->get_source();
        theory_var t = a->get_target();
        numeral  k   = a->get_offset();
        TRACE(assign_profile, tout << "#" << get_enode(s)->get_owner_id() << " #" << get_enode(t)->get_owner_id() << " " << k << "\n";);
        if (l.sign()) {
            k.neg();
            k -= get_epsilon(s);
            add_edge(t, s, k, l);
        }
        else {
            add_edge(s, t, k, l);
        }
        TRACE(ddl_detail, display(tout););
    }

    template<typename Ext>
    void theory_dense_diff_logic<Ext>::new_eq_eh(theory_var v1, theory_var v2) {
        m_arith_eq_adapter.new_eq_eh(v1, v2);
    }

    template<typename Ext>
    bool theory_dense_diff_logic<Ext>::use_diseqs() const {
        return true;
    }

    template<typename Ext>
    void theory_dense_diff_logic<Ext>::new_diseq_eh(theory_var v1, theory_var v2) {
        m_arith_eq_adapter.new_diseq_eh(v1, v2);
    }

    template<typename Ext>
    void theory_dense_diff_logic<Ext>::conflict_resolution_eh(app * atom, bool_var v) {
        // do nothing
    }
    
    template<typename Ext>
    void theory_dense_diff_logic<Ext>::push_scope_eh() {
        theory::push_scope_eh();
        m_scopes.push_back(scope());
        scope & s          = m_scopes.back();
        s.m_atoms_lim      = m_atoms.size();
        s.m_edges_lim      = m_edges.size();
        s.m_cell_trail_lim = m_cell_trail.size();
    }

    template<typename Ext>
    void theory_dense_diff_logic<Ext>::pop_scope_eh(unsigned num_scopes) {
        unsigned lvl     = m_scopes.size();
        SASSERT(num_scopes <= lvl);
        unsigned new_lvl = lvl - num_scopes;
        scope & s        = m_scopes[new_lvl];
        restore_cells(s.m_cell_trail_lim);
        m_edges.shrink(s.m_edges_lim);
        del_atoms(s.m_atoms_lim);
        del_vars(get_old_num_vars(num_scopes));
        m_scopes.shrink(new_lvl);
        theory::pop_scope_eh(num_scopes);
    }

    template<typename Ext>
    void theory_dense_diff_logic<Ext>::restore_cells(unsigned old_size) {
        unsigned sz = m_cell_trail.size();
        unsigned i  = sz;
        while (i > old_size) {
            i--;
            cell_trail & t = m_cell_trail[i];
            cell & c       = m_matrix[t.m_source][t.m_target];
            c.m_edge_id    = t.m_old_edge_id;
            c.m_distance   = t.m_old_distance;
        }
        m_cell_trail.shrink(old_size);
    }

    template<typename Ext>
    void theory_dense_diff_logic<Ext>::del_atoms(unsigned old_size) {
        typename atoms::iterator begin = m_atoms.begin() + old_size;
        typename atoms::iterator it    = m_atoms.end();
        while (it != begin) {
            --it;
            atom * a                      = *it;
            TRACE(del_atoms, tout << "deleting: p" << a->get_bool_var() << "\n";);
            m_bv2atoms[a->get_bool_var()] = 0;
            theory_var s = a->get_source();
            theory_var t = a->get_target();
            TRACE(del_atoms, tout << "m_matrix.size() " << m_matrix.size() << 
                  ", m_matrix[s].size() " << m_matrix[s].size() <<
                  ", m_matrix[t].size(): " << m_matrix[t].size() <<
                  ", t: " << t << ", s: " << s << "\n";);
            SASSERT(m_matrix[s][t].m_occs.back() == a);
            SASSERT(m_matrix[t][s].m_occs.back() == a);
            m_matrix[s][t].m_occs.pop_back();
            m_matrix[t][s].m_occs.pop_back();
            dealloc(a);
        } 
        m_atoms.shrink(old_size);
    }

    template<typename Ext>
    void theory_dense_diff_logic<Ext>::del_vars(unsigned old_num_vars) {
        int num_vars = get_num_vars();
        SASSERT(num_vars >= static_cast<int>(old_num_vars));
        if (num_vars != static_cast<int>(old_num_vars)) {
            m_is_int.shrink(old_num_vars);
            m_f_targets.shrink(old_num_vars);
            m_matrix.shrink(old_num_vars);
            for (auto& cells : m_matrix) {
                cells.shrink(old_num_vars);
            }
        }
    }
        
    template<typename Ext>
    void theory_dense_diff_logic<Ext>::restart_eh() {
        m_arith_eq_adapter.restart_eh();
    }

    template<typename Ext>
    void theory_dense_diff_logic<Ext>::init_search_eh() {
        m_arith_eq_adapter.init_search_eh();
    }

    template<typename Ext>
    final_check_status theory_dense_diff_logic<Ext>::final_check_eh(unsigned level) {
        init_model();
        if (assume_eqs(m_var_value_table))
            return FC_CONTINUE;
        // logical context contains arithmetic expressions that are not
        // in the difference logic fragment.
        if (m_non_diff_logic_exprs)
            return FC_GIVEUP; 
        return FC_DONE;
    }
        
    template<typename Ext>
    bool theory_dense_diff_logic<Ext>::can_propagate() {
        // do nothing
        return false; 
    }

    template<typename Ext>
    void theory_dense_diff_logic<Ext>::propagate() {
        // do nothing
    }

    template<typename Ext>
    void theory_dense_diff_logic<Ext>::flush_eh() {
        // do nothing
    }

    template<typename Ext>
    void theory_dense_diff_logic<Ext>::reset_eh() {
        del_atoms(0);
        m_atoms      .reset();
        m_bv2atoms   .reset();
        m_edges      .reset();
        m_matrix     .reset();
        m_is_int     .reset();
        m_f_targets  .reset();
        m_cell_trail .reset();
        m_scopes     .reset();
        m_non_diff_logic_exprs = false;
        m_edges.push_back(edge());
        theory::reset_eh();
    }
    
    /**
       \brief Store in results the antecedents that justify that the distance between source and target.
    */
    template<typename Ext>
    void theory_dense_diff_logic<Ext>::get_antecedents(theory_var source, theory_var target, literal_vector & result) {
        TRACE(ddl, tout << "get_antecedents, source: #" << get_enode(source)->get_owner_id() << ", target: #" << get_enode(target)->get_owner_id() << "\n";);
        CTRACE(ddl, !is_connected(source, target), display(tout););
        SASSERT(is_connected(source, target));
        svector<var_pair> & todo = m_tmp_pairs;
        todo.reset();

        if (source != target)
            todo.push_back(var_pair(source, target));

        while (!todo.empty()) {
            var_pair & curr = todo.back();
            theory_var s    = curr.first;
            theory_var t    = curr.second;
            todo.pop_back();
        
            SASSERT(is_connected(s, t));
            cell & c        = m_matrix[s][t];
            SASSERT(c.m_edge_id != self_edge_id);
            
            edge & e    = m_edges[c.m_edge_id];
            if (e.m_justification != null_literal)
                result.push_back(e.m_justification);
            
            if (s != e.m_source)
                todo.push_back(var_pair(s, e.m_source));
            if (e.m_target != t)
                todo.push_back(var_pair(e.m_target, t));
        }
    }

    template<typename Ext>
    void theory_dense_diff_logic<Ext>::update_cells() {
        edge_id new_edge_id = m_edges.size() - 1;
        edge & last         = m_edges.back();
        theory_var s        = last.m_source;
        theory_var t        = last.m_target;
        numeral const & k   = last.m_offset;

        // Compute set F of nodes such that:
        // x in F iff
        //    k + d(t, x) < d(s, x)
        
        numeral new_dist;
        row & t_row                = m_matrix[t];
        typename row::iterator it           = t_row.begin();
        typename row::iterator end          = t_row.end();
        typename f_targets::iterator fbegin = m_f_targets.begin();
        typename f_targets::iterator target = fbegin;
        for (theory_var x = 0; it != end; ++it, ++x) {
            if (it->m_edge_id != null_edge_id && x != s) {
                new_dist    = k;
                new_dist   += it->m_distance;
                cell & s_x  = m_matrix[s][x];
                TRACE(ddl, 
                      tout << "s: #" << get_enode(s)->get_owner_id() << " x: #" << get_enode(x)->get_owner_id() << " new_dist: " << new_dist << "\n";
                      tout << "already has edge: " << s_x.m_edge_id << "  old dist: " << s_x.m_distance << "\n";);
                if (s_x.m_edge_id == null_edge_id || new_dist < s_x.m_distance) {
                    target->m_target       = x;
                    target->m_new_distance = new_dist;
                    ++target;
                }
            }
        }
        
        typename f_targets::iterator fend = target;
        
        // For each node y such that y --> s, and for each node x in F,
        // check whether d(y, s) + new_dist(x) < d(y, x).
        typename matrix::iterator it2    = m_matrix.begin();
        typename matrix::iterator end2   = m_matrix.end();
        for (theory_var y = 0; it2 != end2; ++it2, ++y) {
            if (y != t) {
                row  & r = *it2;
                cell & c = r[s];
                if (c.m_edge_id != null_edge_id) {
                    numeral const & d_y_s = c.m_distance;
                    target = fbegin;
                    for (; target != fend; ++target) {
                        theory_var x = target->m_target;
                        if (x != y) {
                            new_dist  = d_y_s;
                            new_dist += target->m_new_distance;
                            cell & y_x = m_matrix[y][x];
                            if (y_x.m_edge_id == null_edge_id || new_dist < y_x.m_distance) {
                                m_cell_trail.push_back(cell_trail(y, x, y_x.m_edge_id, y_x.m_distance));
                                y_x.m_edge_id  = new_edge_id;
                                y_x.m_distance = new_dist;
                                if (!y_x.m_occs.empty()) {
                                    propagate_using_cell(y, x);
                                }
                            }
                        }
                    }
                }
            }
        }
        CASSERT("ddl", check_matrix());
    }

    template<typename Ext>
    void theory_dense_diff_logic<Ext>::assign_literal(literal l, theory_var source, theory_var target) {
        literal_vector & antecedents = m_tmp_literals;
        antecedents.reset();
        get_antecedents(source, target, antecedents);
        ctx.assign(l, ctx.mk_justification(theory_propagation_justification(get_id(), ctx, antecedents.size(), antecedents.data(), l)));
    }
    
    template<typename Ext>
    void theory_dense_diff_logic<Ext>::propagate_using_cell(theory_var source, theory_var target) {
        cell & c = m_matrix[source][target];
        SASSERT(c.m_edge_id != null_edge_id);
        numeral neg_dist = c.m_distance;
        neg_dist.neg();
        typename atoms::const_iterator it  = c.m_occs.begin();
        typename atoms::const_iterator end = c.m_occs.end();
        for (; it != end; ++it) {
            atom * a = *it;
            if (ctx.get_assignment(a->get_bool_var()) == l_undef) {
                if (a->get_source() == source) {
                    SASSERT(a->get_target() == target);
                    if (c.m_distance <= a->get_offset()) {
                        m_stats.m_num_propagations++;
                        TRACE(ddl, tout << "asserting atom to true: "; display_atom(tout, a);
                              tout << "distance(#" << get_enode(source)->get_owner_id() << ", #" << get_enode(target)->get_owner_id()
                              << "): " << c.m_distance << "\n";);
                        assign_literal(literal(a->get_bool_var(), false), source, target);
                    }
                }
                else {
                    SASSERT(a->get_source() == target);
                    SASSERT(a->get_target() == source);
                    if (neg_dist > a->get_offset()) {
                        m_stats.m_num_propagations++;
                        TRACE(ddl, tout << "asserting atom to true: "; display_atom(tout, a);
                              tout << "distance(#" << get_enode(source)->get_owner_id() << ", #" << get_enode(target)->get_owner_id()
                              << "): " << c.m_distance << "\n";);
                        assign_literal(literal(a->get_bool_var(), true), source, target);
                    }
                }
            }
        }
    }

    template<typename Ext>
    inline void theory_dense_diff_logic<Ext>::add_edge(theory_var source, theory_var target, numeral const & offset, literal l) {
        TRACE(ddl, tout << "trying adding edge: #" << get_enode(source)->get_owner_id() << " -- " << offset << " --> #" << get_enode(target)->get_owner_id() << "\n";);
        cell & c_inv = m_matrix[target][source];
        if (c_inv.m_edge_id != null_edge_id && - c_inv.m_distance > offset) {
            // conflict detected.
            TRACE(ddl, tout << "conflict detected: #" << get_enode(source)->get_owner_id() << " #" << get_enode(target)->get_owner_id() <<
                  " offset: " << offset << ", c_inv.m_edge_id: " << c_inv.m_edge_id << ", c_inv.m_distance: " << c_inv.m_distance << "\n";);
            literal_vector & antecedents = m_tmp_literals;
            antecedents.reset();
            get_antecedents(target, source, antecedents);
            if (l != null_literal)
                antecedents.push_back(l);
            ctx.set_conflict(ctx.mk_justification(theory_conflict_justification(get_id(), ctx, antecedents.size(), antecedents.data())));

            return;
        }
        
        cell & c = m_matrix[source][target];
        if (c.m_edge_id == null_edge_id || offset < c.m_distance) {
            TRACE(ddl, tout << "adding edge: #" << get_enode(source)->get_owner_id() << " -- " << offset << " --> #" << get_enode(target)->get_owner_id() << "\n";);
            m_edges.push_back(edge(source, target, offset, l));
            update_cells();
        }
    }


#ifdef Z3DEBUG
    template<typename Ext>
    bool theory_dense_diff_logic<Ext>::check_vector_sizes() const {
        SASSERT(m_matrix.size() == m_f_targets.size());
        SASSERT(m_is_int.size() == m_matrix.size());
        typename matrix::const_iterator it  = m_matrix.begin();
        typename matrix::const_iterator end = m_matrix.end();
        for (; it != end; ++it) {
            SASSERT(it->size() == m_matrix.size());
        }
        return true;
    }

    template<typename Ext>
    bool theory_dense_diff_logic<Ext>::check_matrix() const {
        int sz = static_cast<int>(m_matrix.size());
        for (theory_var i = 0; i < sz; ++i) {
            for (theory_var j = 0; j < sz; ++j) {
                cell const & c = m_matrix[i][j];
                if (c.m_edge_id == self_edge_id) {
                    SASSERT(i == j);
                    SASSERT(c.m_distance.is_zero());
                }
                else if (c.m_edge_id != null_edge_id) {
                    edge const & e = m_edges[c.m_edge_id];
                    theory_var s = e.m_source;
                    theory_var t = e.m_target;
                    numeral k  = get_distance(i, s);
                    k         += e.m_offset;
                    k         += get_distance(t, j);
                    if (c.m_distance != k) {
                        CTRACE(ddl, c.m_distance != k, tout << "i: " << i << " j: " << j << " k: " << k << " c.m_distance: " << c.m_distance << "\n";
                              display(tout););
                        SASSERT(c.m_distance == k);
                    }
                }
            }
        }
        return true;
    }
#endif        

    template<typename Ext>
    void theory_dense_diff_logic<Ext>::display(std::ostream & out) const {
        out << "Theory dense difference logic:\n";
        display_var2enode(out);
        typename matrix::const_iterator it1  = m_matrix.begin();
        typename matrix::const_iterator end1 = m_matrix.end();
        for (int v1 = 0; it1 != end1; ++it1, ++v1) {
            typename row::const_iterator it2  = it1->begin();
            typename row::const_iterator end2 = it1->end();
            for (int v2 = 0; it2 != end2; ++it2, ++v2) {
                if (it2->m_edge_id != null_edge_id && it2->m_edge_id != self_edge_id) {
                    out << "#";
                    out.width(5);
                    out << std::left << get_enode(v1)->get_owner_id() << " -- ";
                    out.width(10);
                    out << std::left << it2->m_distance << " : id";
                    out.width(5);
                    out << std::left << it2->m_edge_id << " --> #";
                    out << get_enode(v2)->get_owner_id() << "\n";
                }
            }
        }
        out << "atoms:\n";
        typename atoms::const_iterator it2  = m_atoms.begin();
        typename atoms::const_iterator end2 = m_atoms.end();
        for (;it2 != end2; ++it2) {
            atom * a = *it2;
            display_atom(out, a);
        }
    }

    template<typename Ext>
    void theory_dense_diff_logic<Ext>::display_atom(std::ostream & out, atom * a) const {
        out << "#";
        out.width(5);
        out << std::left << get_enode(a->get_target())->get_owner_id() << " - #";
        out.width(5);
        out << std::left << get_enode(a->get_source())->get_owner_id() << " <= ";
        out.width(10);
        out << std::left << a->get_offset() << "        assignment: " << ctx.get_assignment(a->get_bool_var()) << "\n";
    }

    template<typename Ext>
    void theory_dense_diff_logic<Ext>::collect_statistics(::statistics & st) const {
        st.update("dd assertions", m_stats.m_num_assertions);
        st.update("dd propagations", m_stats.m_num_propagations);
        m_arith_eq_adapter.collect_statistics(st);
    }

    /**
       \brief Build a model for doing model-based theory combination.
    */
    template<typename Ext>
    void theory_dense_diff_logic<Ext>::init_model() {
        int num_vars = get_num_vars();
        m_assignment.reset();
        m_assignment.resize(num_vars);
        for (int i = 0; i < num_vars; ++i) {
            row & r     = m_matrix[i];
            numeral & d = m_assignment[i];
            for (int j = 0; j < num_vars; ++j) {
                if (i != j) {
                    cell & c = r[j];
                    if (c.m_edge_id != null_edge_id && c.m_distance < d) {
                        d = c.m_distance;
                    }
                }
            }
        }
        for (int i = 0; i < num_vars; ++i)
            m_assignment[i].neg();
        TRACE(ddl_model, 
              tout << "ddl model\n";
              for (theory_var v = 0; v < num_vars; ++v) {
                  tout << "#" << mk_pp(get_enode(v)->get_expr(), m) << " = " << m_assignment[v] << "\n";
              });
    }

    /**
       The arithmetic module uses infinitesimals. So,
       an inf_numeral (n,k) represents  n + k*epsilon
       where epsilon is a very small number. 
       In order to generate a model, we need to compute
       a value for epsilon in a way all inequalities remain
       satisfied.

       assume we have the inequality
       x - y <= (n_c, k_c)

       where the interpretation for x and y is:
       (n_x, k_x) and (n_y, k_y).

       So,
       
       (n_x, k_x) <= (n_y + n_c, k_y + k_c)
       

       The only interesting case is n_x < n_y + n_c and k_x > k_y + k_c.
       Using the definition of infinitesimal numbers
       we have:
       
       n_x + k_x * epsilon <= n_y + n_c + (k_y + k_c) * epsilon

       Therefore:
       
       epsilon <= (n_y + n_c - n_x) / 2*(k_x - k_y - k_c)

       Division by 2 is introduced to avoid introducing new equalities.
       Consider when n_x = n_y = k_x = k_y = 0, n_c = 4, k_c = -1,
       then n_c / -k_c = n_c, hence setting epsilon to n_c introduces a value
       that is not smaller than the smallest integral difference.
    */
    template<typename Ext>
    void theory_dense_diff_logic<Ext>::compute_epsilon() {
        m_epsilon = rational(1, 2);
        TRACE(ddl, display(tout););
        typename edges::const_iterator it  = m_edges.begin();
        typename edges::const_iterator end = m_edges.end();
        // first edge is null
        SASSERT(it->m_target == null_theory_var);
        SASSERT(it->m_source == null_theory_var);
        ++it; 
        for (; it != end; ++it) {
            edge const & e = *it;
            rational n_x = m_assignment[e.m_target].get_rational().to_rational();
            rational k_x = m_assignment[e.m_target].get_infinitesimal().to_rational();
            rational n_y = m_assignment[e.m_source].get_rational().to_rational();
            rational k_y = m_assignment[e.m_source].get_infinitesimal().to_rational();
            rational n_c = e.m_offset.get_rational().to_rational();
            rational k_c = e.m_offset.get_infinitesimal().to_rational();
            TRACE(ddl,
                  tout << e.m_source << " - " << e.m_target << " <= " << e.m_offset << "\n";
                  tout << "(n_x,k_x): " << n_x << ", " << k_x << 
                  ", (n_y,k_y): " << n_y << ", " << k_y << 
                  ", (n_c,k_c): " << n_c << ", " << k_c << "\n";);
            if (n_x < n_y + n_c && k_x > k_y + k_c) {
                rational new_epsilon = (n_y + n_c - n_x) / (2*(k_x - k_y - k_c));
                if (new_epsilon < m_epsilon) {
                    TRACE(ddl, tout << "new epsilon: " << new_epsilon << "\n";);
                    m_epsilon = new_epsilon;
                }
            }
        }
    }

    template<typename Ext>
    void theory_dense_diff_logic<Ext>::fix_zero() {
        int num_vars = get_num_vars();
        for (int v = 0; v < num_vars && v < (int)m_assignment.size(); ++v) {
            enode * n = get_enode(v);
            if (m_autil.is_zero(n->get_expr()) && !m_assignment[v].is_zero()) {
                numeral val = m_assignment[v];
                sort * s = n->get_expr()->get_sort();
                // adjust the value of all variables that have the same sort.
                for (int v2 = 0; v2 < num_vars; ++v2) {
                    enode * n2 = get_enode(v2);
                    if (n2->get_expr()->get_sort() == s) {
                        m_assignment[v2] -= val;
                    }
                }
                SASSERT(m_assignment[v].is_zero());
            }
        }
        TRACE(ddl_model, 
              tout << "ddl model\n";
              for (theory_var v = 0; v < num_vars; ++v) {
                  tout << "#" << mk_pp(get_enode(v)->get_expr(), m) << " = " << m_assignment[v] << "\n";
              });
    }

    template<typename Ext>
    void theory_dense_diff_logic<Ext>::init_model(model_generator & mg) {
        m_factory = alloc(arith_factory, m);
        mg.register_factory(m_factory);
        if (!m_assignment.empty()) {
            fix_zero();
            compute_epsilon();
        }
    }

    template<typename Ext>
    model_value_proc * theory_dense_diff_logic<Ext>::mk_value(enode * n, model_generator & mg) {
        theory_var v = n->get_th_var(get_id());
        SASSERT(v != null_theory_var);
        if (v >= (int)m_assignment.size())
            return alloc(expr_wrapper_proc, m_factory->mk_num_value(rational::zero(), is_int(v)));
        numeral const & val = m_assignment[v];
        rational num = val.get_rational().to_rational() +  m_epsilon *  val.get_infinitesimal().to_rational();
        return alloc(expr_wrapper_proc, m_factory->mk_num_value(num, is_int(v)));
    }

    // TBD: code is common to both sparse and dense difference logic solvers.
    template<typename Ext>
    bool theory_dense_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_autil.is_numeral(n, r)) {
            q += r;
        }
        else if (m_autil.is_add(n)) {
            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_autil.is_mul(n, x, y) && m_autil.is_numeral(x, r)) {
            return internalize_objective(y, m*r, q, objective);
        }
        else if (m_autil.is_mul(n, y, x) && m_autil.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_autil.get_family_id()) {
            return false;
        }
        else {
            enode * e = nullptr;
            theory_var v = 0;
            if (ctx.e_internalized(n)) {
                e = ctx.get_enode(to_app(n));                
            }
            else {
                ctx.internalize(n, false);
                e = ctx.get_enode(n);
            }            
            v = e->get_th_var(get_id());
            if (v == null_theory_var) {
                v = mk_var(e);                
            }
            bool found = false;
            for (auto& p : objective) {
                if (p.first == v) {
                    p.second += m;
                    found = true;
                }
            }
            if (!found) {
                objective.push_back(std::make_pair(v, m));
            }
        }
        return true;
    }

    template<typename Ext>
    inf_eps_rational<inf_rational> theory_dense_diff_logic<Ext>::value(theory_var v) {
        objective_term const& objective = m_objectives[v];   
        inf_eps r = inf_eps(m_objective_consts[v]);
        for (unsigned i = 0; i < objective.size(); ++i) {
            numeral n = m_assignment[v];
            rational r1 = n.get_rational().to_rational();
            rational r2 = n.get_infinitesimal().to_rational();
            r += objective[i].second * inf_eps(rational(0), inf_rational(r1, r2));
        }
        return r;
    }

    template<typename Ext>
    inf_eps_rational<inf_rational> theory_dense_diff_logic<Ext>::maximize(theory_var v, expr_ref& blocker, bool& has_shared) {
        typedef simplex::simplex<simplex::mpq_ext> Simplex;
        Simplex S(m.limit());
        objective_term const& objective = m_objectives[v];
        has_shared = false;
        
        IF_VERBOSE(4,
                   for (auto const& o : objective) {
                       verbose_stream() << o.second << " * v" << o.first << " ";
                   }
                   verbose_stream() << " + " << m_objective_consts[v] << "\n";);

        unsynch_mpq_manager mgr;
        unsynch_mpq_inf_manager inf_mgr;
        unsigned num_nodes = get_num_vars();
        unsigned num_edges = m_edges.size();
        S.ensure_var(num_nodes + num_edges + m_objectives.size());
        for (unsigned i = 0; i < num_nodes; ++i) {
            numeral const& a = m_assignment[i];
            rational fin = a.get_rational().to_rational();
            rational inf = a.get_infinitesimal().to_rational();
            mpq_inf q(mgr.dup(fin.to_mpq()), mgr.dup(inf.to_mpq()));
            S.set_value(i, q);
            inf_mgr.del(q);
        }
        for (unsigned i = 0; i < num_nodes; ++i) {
            enode * n = get_enode(i);
            if (m_autil.is_zero(n->get_expr())) {
                S.set_lower(i, mpq_inf(mpq(0), mpq(0)));
                S.set_upper(i, mpq_inf(mpq(0), mpq(0)));
                break;
            }
        }
        svector<unsigned> vars;
        scoped_mpq_vector coeffs(mgr);
        coeffs.push_back(mpq(1));
        coeffs.push_back(mpq(-1));
        coeffs.push_back(mpq(-1));
        vars.resize(3);
        for (unsigned i = 0; i < num_edges; ++i) {
            edge const& e = m_edges[i];
            if (e.m_source == null_theory_var || e.m_target == null_theory_var) {
                continue;
            }
            unsigned base_var = num_nodes + i;
            vars[0] = e.m_target;
            vars[1] = e.m_source;
            vars[2] = base_var;
            S.add_row(base_var, 3, vars.data(), coeffs.data());
            // t - s <= w
            // =>
            // t - s - b = 0, b <= w
            numeral const& w = e.m_offset;
            rational fin = w.get_rational().to_rational();
            rational inf = w.get_infinitesimal().to_rational();
            mpq_inf q(mgr.dup(fin.to_mpq()), mgr.dup(inf.to_mpq()));
            S.set_upper(base_var, q);            
            inf_mgr.del(q);
        }
        unsigned w = num_nodes + num_edges + v;

        // add objective function as row.
        coeffs.reset();
        vars.reset();
        for (auto const& o : objective) {
            coeffs.push_back(o.second.to_mpq());
            vars.push_back(o.first);
        }
        coeffs.push_back(mpq(1));
        vars.push_back(w);
        Simplex::row row = S.add_row(w, vars.size(), vars.data(), coeffs.data());
        
        TRACE(opt, S.display(tout); display(tout););
        
        // optimize    
        lbool is_sat = S.make_feasible();
        if (is_sat == l_undef) {
            blocker = m.mk_false();
            return inf_eps::infinity();        
        }
        TRACE(opt, S.display(tout); );    
        SASSERT(is_sat != l_false);
        lbool is_fin = S.minimize(w);

        ensure_rational_solution(S);
        
        switch (is_fin) {
        case l_true: {
            simplex::mpq_ext::eps_numeral const& val = S.get_value(w);
            inf_rational r(-rational(val.first), -rational(val.second));
            TRACE(opt, tout << r << " " << "\n"; 
                  S.display_row(tout, row, true););
            Simplex::row_iterator it = S.row_begin(row), end = S.row_end(row);
            expr_ref_vector& core = m_objective_assignments[v];
            expr_ref tmp(m);
            core.reset();
            for (; it != end; ++it) {
                unsigned v = it->var();
                if (num_nodes <= v && v < num_nodes + num_edges) {
                    unsigned edge_id = v - num_nodes;
                    literal lit = m_edges[edge_id].m_justification;
                    if (lit != null_literal) {
                        ctx.literal2expr(lit, tmp);
                        core.push_back(tmp);
                    }
                }
                TRACE(opt, tout << core << "\n";);
            }
            for (unsigned i = 0; i < num_nodes; ++i) {
                mpq_inf const& val = S.get_value(i);
                rational q(val.first);
                SASSERT(rational(val.second).is_zero());
                numeral  a(q);
                m_assignment[i] = a;
            }
            inf_eps result(rational(0), r);
            blocker = mk_gt(v, result);
            IF_VERBOSE(10, verbose_stream() << blocker << "\n";);
            r += m_objective_consts[v];
            return inf_eps(rational(0), r);
        }
        default:
            TRACE(opt, tout << "unbounded\n"; );        
            blocker = m.mk_false();
            return inf_eps::infinity();        
        }
    }

    template<typename Ext>
    theory_var theory_dense_diff_logic<Ext>::add_objective(app* term) {
        TRACE(opt, tout << mk_pp(term, m) << "\n";);
        objective_term objective;
        theory_var result = m_objectives.size();
        rational q(1), r(0);
        expr_ref_vector vr(m);
        if (!is_linear(m, term)) {
            result = null_theory_var;
        }
        else if (internalize_objective(term, q, r, objective)) {
            m_objectives.push_back(objective);
            m_objective_consts.push_back(r);
            m_objective_assignments.push_back(vr);
        }
        else {
            result = null_theory_var;
        }
        return result;
    }

    template<typename Ext>
    expr_ref theory_dense_diff_logic<Ext>::mk_gt(theory_var v, inf_eps const& val) {
        return mk_ineq(v, val, true);
    }

    template<typename Ext>
    expr_ref theory_dense_diff_logic<Ext>::mk_ge(
        generic_model_converter& fm, theory_var v, inf_eps const& val) {
        return mk_ineq(v, val, false);
    }

    template<typename Ext>
    expr_ref theory_dense_diff_logic<Ext>::mk_ineq(theory_var v, inf_eps const& val, bool is_strict) {
        objective_term const& t = m_objectives[v];
        expr_ref e(m), f(m), f2(m);
        TRACE(opt, tout << "mk_ineq " << v << " " << val << "\n";);
        if (t.size() == 1 && t[0].second.is_one()) {
            f = get_enode(t[0].first)->get_expr();
        }
        else if (t.size() == 1 && t[0].second.is_minus_one()) {
            f = m_autil.mk_uminus(get_enode(t[0].first)->get_expr());
        }
        else if (t.size() == 2 && t[0].second.is_one() && t[1].second.is_minus_one()) {
            f = get_enode(t[0].first)->get_expr();
            f2 = get_enode(t[1].first)->get_expr();
            f = m_autil.mk_sub(f, f2); 
        }
        else if (t.size() == 2 && t[1].second.is_one() && t[0].second.is_minus_one()) {
            f = get_enode(t[1].first)->get_expr();
            f2 = get_enode(t[0].first)->get_expr();
            f = m_autil.mk_sub(f, f2);
        }
        else {
            // 
            expr_ref_vector const& core = m_objective_assignments[v];
            f = m.mk_and(core);
            if (is_strict) {
                f = m.mk_not(f);
            }
            TRACE(arith, tout << "block: " << f << "\n";);
            return f;
        }
        
        e = m_autil.mk_numeral(val.get_rational(), f->get_sort());
        
        if (val.get_infinitesimal().is_neg()) {
            if (is_strict) {
                f = m_autil.mk_ge(f, e);
            }
            else {
                expr_ref_vector const& core = m_objective_assignments[v];
                f = m.mk_and(core);
            }
        }
        else {
            if (is_strict) {
                f = m_autil.mk_gt(f, e);
            }
            else {
                f = m_autil.mk_ge(f, e);
            }
        }
        return f;
    }

};