z3-src 503.2.0

/*++
Copyright (c) 2013 Microsoft Corporation

Module Name:

    optsmt.cpp

Abstract:
   
    Objective optimization method.

Author:

    Anh-Dung Phan (t-anphan) 2013-10-16

Notes:

    Suppose we obtain solution t1 = k1, ..., tn = kn-epsilon
    Assert:
      t1 > k1 \/ t2 > k2 \/ ... \/ tn >= kn
    If this solution is satisfiable, then for each t_i, maximize the 
    assignment and assert the new frontier.
    Claim: we don't necessarily have to freeze assignments of 
    t_i when optimizing assignment for t_j
    because the state will always satisfy the disjunction.
    If one of the k_i is unbounded, then omit a disjunction for it.        
    

--*/

#include <typeinfo>
#include <cstring>
#include <optional>
#include "util/common_msgs.h"
#include "opt/optsmt.h"
#include "opt/opt_geometric.h"
#include "opt/opt_search.h"
#include "opt/opt_nlsat.h"
#include "opt/opt_solver.h"
#include "opt/opt_context.h"
#include "ast/arith_decl_plugin.h"
#include "smt/theory_arith.h"
#include "ast/ast_pp.h"
#include "ast/ast_util.h"
#include "model/model_pp.h"
#include "ast/rewriter/th_rewriter.h"
#include "opt/opt_params.hpp"

namespace opt {

    namespace {
        // Track pushes used to try temporary bounds such as obj >= best + step.
        // These bounds may be infeasible, so remove them before nlsat/bisection
        // and on exit.
        class scoped_pushes {
            solver& m_solver;
            unsigned m_size = 0;
        public:
            explicit scoped_pushes(solver& s): m_solver(s) {}
            ~scoped_pushes() { reset(); }
            scoped_pushes(scoped_pushes const&) = delete;
            scoped_pushes& operator=(scoped_pushes const&) = delete;

            unsigned size() const { return m_size; }

            void push() {
                m_solver.push();
                ++m_size;
            }

            void pop() {
                SASSERT(m_size > 0);
                m_solver.pop(1);
                --m_size;
            }

            void reset() {
                if (m_size > 0)
                    m_solver.pop(m_size);
                m_size = 0;
            }
        };
    }

    void optsmt::set_max(vector<objective_value>& dst, vector<inf_eps> const& src, expr_ref_vector& fmls) {
        for (unsigned i = 0; i < src.size(); ++i) {
            if (src[i] >= dst[i].rational_bound()) {
                dst[i] = src[i];
                m_models.set(i, m_s->get_model_idx(i));
                m_s->get_labels(m_labels);
                m_lower_fmls[i] = fmls.get(i);
                if (dst[i].rational_bound().is_pos() && !dst[i].is_finite()) { // review: likely done already.
                    m_lower_fmls[i] = m.mk_false();
                    fmls[i] = m.mk_false();
                }
            }
            else if (src[i] < dst[i].rational_bound() && !m.is_true(m_lower_fmls.get(i))) {
                fmls[i] = m_lower_fmls.get(i);                
            }
        }
    }

    /*
        Enumerate locally optimal assignments until fixedpoint.
    */
    lbool optsmt::geometric_opt() {
        lbool is_sat = l_true;

        expr_ref bound(m), last_bound(m);

        vector<inf_eps> previous_lower;
        for (auto const& value : m_lower)
            previous_lower.push_back(value.rational_bound());
        geometric_step step;
        scoped_pushes scopes(*m_s);
        unsigned delta_index = 0;    // index of objective to speed up.
        bool has_bound = false;      // is the current objective bounded by a constraint.

        while (m.inc()) {
            is_sat = m_s->check_sat(0, nullptr);
            if (is_sat == l_true) { 
                bound = update_lower();
                if (!m.is_true(bound))
                    has_bound = true;
                step.update(can_increment_delta(previous_lower, delta_index));
                if (step.value() > rational::one()) {
                    scopes.push();
                    // only try to improve delta_index. 
                    bound = m_s->mk_ge(delta_index, lower(delta_index) + inf_eps(step.value()));
                }
                TRACE(opt, tout << mk_pp(m_objs.get(delta_index), m) << " index: " << delta_index
                      << " delta: " << step.value() << " bound: " << bound
                      << " " << lower(delta_index) << " " << upper(delta_index) << "\n");
                if (bound == last_bound) {
                    is_sat = l_false;
                    if ((!has_bound || !m_lower[delta_index].is_finite()) && !m_upper[delta_index].is_finite())
                        m_lower[delta_index] = m_upper[delta_index];
                }
                else {
                    m_s->assert_expr(bound);                
                    last_bound = bound;                    
                    continue;
                }
            }
            if (is_sat == l_false && step.value() > rational::one()) {
                step.reset();
                scopes.pop();
                last_bound = nullptr;
            }
            else if (is_sat == l_false) {
                // we are done with this delta_index.
                m_upper[delta_index] = m_lower[delta_index];
                scopes.reset();
                last_bound = nullptr;
                bool all_tight = true;
                for (unsigned i = 0; i < m_lower.size(); ++i) {
                    all_tight &= lower(i) == upper(i);
                }
                if (all_tight || delta_index + 1 == m_lower.size())
                    break;
                step.reset();
                ++delta_index;
                has_bound = false;
            }
            else {
                scopes.reset();
                break;
            }
        }
        
        if (!m.inc() || is_sat == l_undef) {
            return l_undef;
        }
        
        return l_true;        
    }

    bool optsmt::is_unbounded(unsigned obj_index, bool is_maximize) {
        if (is_maximize) {
            return !m_upper[obj_index].is_finite();
        }
        else {
            return !m_lower[obj_index].is_finite();
        }
    }

    lbool optsmt::geometric_search(unsigned obj_index, bool is_maximize, unsigned unbounded_check_rounds) {
        TRACE(opt, tout << "index: " << obj_index << " is-max: " << is_maximize << "\n";);
        arith_util arith(m);
        bool is_int = arith.is_int(m_objs.get(obj_index));
        lbool is_sat = l_true;
        expr_ref bound(m), last_bound(m);

        geometric_step step;
        scoped_pushes scopes(*m_s);
        objective_search search(upper(obj_index));
        inf_eps last_objective = inf_eps(rational(-1), inf_rational(0));
        inf_eps const infty(rational(1), inf_rational(0));
        bool last_bound_valid = true;
        unsigned climb_rounds = 0;
        // initial budget: same constant as the bounded nlsat run in
        // nla_core::bounded_nlsat; the known unbounded proofs need 4k-23k
        unsigned unbounded_check_rlimit = 100000;

        while (m.inc()) {
            is_sat = m_s->check_sat(0, nullptr);
            TRACE(opt, tout << "check " << is_sat << "\n";
                  tout << "last bound: " << last_bound << " bound " << bound << "\n";
                  tout << "lower: " << lower(obj_index) << "\n";
                  tout << "upper: " << upper(obj_index) << "\n";
                  if (is_sat == l_true) m_s->display(tout);
                  );
            if (is_sat == l_true) {                
                // Once a finite cap is known, leave refinement to the normal
                // search instead of perturbing it with more certification probes.
                bool probe_bound = !search.upper().is_finite();
                auto result = m_s->maximize_objective(obj_index, bound, m_dual_bounds && !is_int, probe_bound);
                last_bound_valid = result.bound_valid;
                search.update(result);
                m_upper[obj_index] = search.upper();
                m_s->get_model(m_model);
                SASSERT(m_model);
                inf_eps obj = m_s->saved_objective_value(obj_index);
                TRACE(opt, tout << "saved objective: " << obj << "\n";);
                update_lower_lex(obj_index, obj, is_maximize);
                if (m.limit().is_canceled()) {
                    is_sat = l_undef;
                    break;
                }
                step.update(is_int && m_lower[obj_index].is_finite());
                // A real objective may improve forever without refuting an
                // upper bound. After a streak, use prove_unbounded_above to
                // commit +oo with the current model as a witness. Bounded
                // objectives can show the same pattern, so limit the query by
                // deterministic resource count and exponentially back off
                // inconclusive retries.
                if (is_int || search.upper().is_finite())
                    climb_rounds = 0;
                else if (m_optsmt_nlsat && ++climb_rounds >= unbounded_check_rounds) {
                    climb_rounds = 0;
                    unbounded_check_rounds *= 2;
                    if (prove_unbounded_above(obj_index, unbounded_check_rlimit)) {
                        set_best(obj_index, infty, is_maximize);
                        break;
                    }
                    // saturate: an overflow to 0 would make the budget
                    // infinite (reslimit::push treats 0 as no limit)
                    if (unbounded_check_rlimit <= UINT_MAX / 2)
                        unbounded_check_rlimit *= 2;
                }
                // When maximize_objective could not validate its arithmetic
                // hint (result.bound_valid == false), the blocker refers to
                // that unachievable hint and must not be used.  'obj' now holds
                // the value of an actual model, so replace the blocker with a
                // model-derived tightening so the search keeps making progress
                // toward the true optimum instead of terminating prematurely
                // (issue #10028).
                if (!result.bound_valid || step.value() > rational::one() || (obj == last_objective && is_int)) {
                    scopes.push();
                    bound = m_s->mk_ge(obj_index, obj + inf_eps(step.value()));
                    search.enter_region(obj + inf_eps(step.value()));
                }
                last_objective = obj;
                if (bound == last_bound) {
                    // LP didn't produce a new blocker. If the model-based lower bound
                    // is strictly better than what the LP found, use it to push the LP
                    // further. This handles cases where nonlinear constraints, mod,
                    // to_int, prevent the LP from seeing the full feasible region.
                    if (m_lower[obj_index].is_finite() && lower(obj_index) > obj)
                        bound = m_s->mk_ge(obj_index, lower(obj_index));
                    if (bound == last_bound && obj.get_infinitesimal().is_pos()) {
                        // The objective sits infinitesimally above the strict
                        // bound asserted in the previous round: r + k*delta with
                        // k > 0.  Such a value is not a proven optimum, it only
                        // says the arithmetic solver could not move off the
                        // bound it was just given.  Its blocker is built from
                        // the rational part alone, so it collapses onto the
                        // previous blocker and the search would stop here and
                        // report the stalled value as the optimum.  Force a
                        // strictly larger rational step instead; if the step is
                        // infeasible the loop terminates through the l_false
                        // branch below with the best proven bound.
                        scopes.push();
                        bound = m_s->mk_ge(obj_index, obj + inf_eps(step.value()));
                        search.enter_region(obj + inf_eps(step.value()));
                    }
                    if (bound == last_bound)
                        break;
                }
                // A validated blocker can discard better values than the
                // enclosing trial cut without introducing another scope.
                search.exclude_below(std::max(obj, lower(obj_index)));
                m_s->assert_expr(bound);
                last_bound = bound;
            }
            else if (is_sat == l_false && step.value() > rational::one()) {
                step.reset();
                scopes.pop();
                search.leave_region(is_sat);
                m_upper[obj_index] = search.upper();
            }
            else {
                if (is_sat == l_false) {
                    search.refute_region();
                    m_upper[obj_index] = search.upper();
                }
                if (((is_sat == l_false && !last_bound_valid) || is_sat == l_undef) && !is_int && m_lower[obj_index].is_finite()) {
                    // A real-valued objective whose last step was the
                    // model-derived fallback: the failed step only shows
                    // the trial region is empty, it does not make the model value
                    // optimal. Close or narrow the gap by bisection instead of
                    // reporting the lower bound as the optimum.
                    inf_eps hi = search.upper();
                    if (lower(obj_index) < hi) {
                        scopes.reset();
                        search.reset_regions();
                        bool smt_gave_up = is_sat == l_undef;
                        is_sat = refine_real_objective(obj_index, is_maximize, hi, smt_gave_up);
                    }
                }
                break;
            }
        }
        TRACE(opt, tout << is_sat << " " << scopes.size() << "\n";);
        scopes.reset();

        if (is_sat == l_false && !m_model) {
            return l_false;
        }
        
        if (!m.inc() || is_sat == l_undef) {
            return l_undef;
        }

        // A satisfiable objective whose lower bound is still -oo means no
        // model value was ever extracted before the loop stalled. Tightening
        // below would then report -oo as the optimum of a satisfiable
        // objective; report the interval as unknown instead.
        if (!m_lower[obj_index].is_finite() && lower(obj_index).is_neg())
            return l_undef;

        // Exact endpoints share a finite part but retain their outward
        // rational brackets for the legacy search.
        if (!m_lower[obj_index].exact_finite())
            m_upper[obj_index] = m_lower[obj_index];
        return l_true;
    }

    void optsmt::publish_best_model() {
        // Callbacks may add constraints, so save the model and labels first.
        m_best_model = m_model;
        m_s->get_labels(m_labels);
        m_context.set_model(m_model);
    }

    void optsmt::set_best(unsigned idx, inf_eps const& v, bool is_maximize) {
        m_lower[idx] = v;
        m_upper[idx].reset_exact();
        IF_VERBOSE(1, 
                   if (is_maximize) 
                       verbose_stream() << "(optsmt lower bound: " << v << ")\n";
                   else
                       verbose_stream() << "(optsmt upper bound: " << (-v) << ")\n";
                   );
        publish_best_model();
    }

    /**
       \brief Certify that objective idx is unbounded above over the current
       assertions by projecting feasible nlsat regions onto the objective,
       given a budget of rlimit_budget resource-counter ticks; bounds pushed
       by the climb hold in the current model and are harmless. The budget's
       expiry surfaces from the nlsat solver as a solver_exception, not a
       tactic_exception, so any z3_exception here means "no proof within
       budget", never an error to propagate. A cancellation flag (a soft
       timeout, ctrl-c) also lands in the same catch, and popping the budget
       would erase it; it is probed while still set and re-raised so the
       caller's m.inc() loop still observes it. An expired tighter enclosing
       rlimit needs no such care: the counter stays above that ceiling after
       the pop and keeps failing on its own.
    */
    bool optsmt::prove_unbounded_above(unsigned idx, unsigned rlimit_budget) {
        if (!m_lower[idx].is_finite())
            return false;
        expr_ref_vector hard(m);
        for (unsigned i = 0; i < m_s->get_num_assertions(); ++i)
            hard.push_back(m_s->get_assertion(i));
        params_ref p;
        nlsat_opt engine(m, p);
        lbool r = l_undef;
        uint64_t start = m.limit().count();
        uint64_t consumed = 0;
        bool external_cancel = false;
        {
            scoped_rlimit budget(m.limit(), rlimit_budget);
            try {
                r = engine.prove_unbounded(hard, m_objs.get(idx), lower(idx).get_rational());
            }
            catch (z3_exception&) {
                external_cancel = m.limit().get_cancel_msg() == Z3_CANCELED_MSG;
                r = l_undef;
            }
            // read before the pop of the budget truncates the counter
            consumed = m.limit().count() - start;
        }
        if (external_cancel)
            m.limit().cancel();
        IF_VERBOSE(1, if (r == l_true) verbose_stream() << "(optsmt nlsat unbounded above)\n");
        IF_VERBOSE(3, verbose_stream() << "(optsmt unbounded-check :result " << r
                   << " :rlimit-consumed " << consumed
                   << " :budget " << rlimit_budget << ")\n");
        TRACE(opt, tout << "prove_unbounded_above: " << r << " budget: " << rlimit_budget << "\n";);
        return r == l_true;
    }

    lbool optsmt::refine_real_objective(unsigned idx, bool is_maximize, inf_eps const& hi, bool smt_gave_up) {
        if (m_optsmt_nlsat) {
            nlsat_outcome outcome;
            try {
                outcome = nlsat_cells(idx, is_maximize, hi);
            }
            catch (default_exception const& ex) {
                if (!m.limit().is_canceled() ||
                    (std::strcmp(ex.what(), Z3_CANCELED_MSG) != 0 && std::strcmp(ex.what(), Z3_MAX_RESOURCE_MSG) != 0))
                    throw;
                // Algebraic operations use "canceled" for resource exhaustion too.
                // Preserve the proven interval, but do not swallow unrelated errors.
                m_s->set_reason_unknown(m.limit().get_cancel_msg());
                return l_undef;
            }
            switch (outcome) {
            case nlsat_outcome::certified:
                return l_true;
            case nlsat_outcome::incomplete:
                return l_undef;
            case nlsat_outcome::no_result:
                break;
            }
        }
        if (!hi.is_finite() || smt_gave_up)
            return l_undef;
        if (lower(idx) < hi)
            return bisect(idx, is_maximize, hi);
        return l_true;
    }

    /**
       \brief Exact optimization over nlsat cells:
       maximize the objective over the hard constraints within
       [m_lower[idx], hi] with nlsat_opt. On success the optimum may be an
       algebraic number: both objective_value endpoints carry its exact
       finite part and outward rational bounds, with -epsilon for an open limit.
       Returns certified for an attained optimum, an open limit, or unboundedness;
       incomplete when an accepted model leaves an unresolved interval; and
       no_result when no result is accepted (unsupported fragment, no model,
       or a lower bracket below the incumbent). Only no_result permits the
       caller to try another fallback instead of retaining the nlsat result.
    */
    optsmt::nlsat_outcome optsmt::nlsat_cells(unsigned idx, bool is_maximize, inf_eps const& hi) {
        if (!m_lower[idx].is_finite())
            return nlsat_outcome::no_result;
        expr_ref_vector hard(m);
        for (unsigned i = 0; i < m_s->get_num_assertions(); ++i)
            hard.push_back(m_s->get_assertion(i));
        params_ref p;
        nlsat_opt engine(m, p);
        nlsat_opt::result res(m);
        IF_VERBOSE(2, verbose_stream() << "(optsmt nlsat [" << lower(idx) << ", " << hi << "])\n");
        std::optional<rational> rational_upper;
        if (hi.is_finite())
            rational_upper = hi.get_rational();
        lbool r = engine.maximize(hard, m_objs.get(idx), lower(idx).get_rational(), rational_upper,
                                 m_bisect_rounds, res, m_nlsat_supremum_rlimit);
        TRACE(opt, tout << "nlsat cells: " << r << " rounds " << res.m_rounds << " value " << res.m_value << "\n";);
        if (!res.m_model)
            return nlsat_outcome::no_result;
        if (res.m_unbounded) {
            m_model = res.m_model;
            set_best(idx, inf_eps(rational(1), inf_rational(0)), is_maximize);
            m_upper[idx] = m_lower[idx];
            return nlsat_outcome::certified;
        }
        if (res.m_open) {
            // The real model is a feasible witness, not an assignment at the
            // limit. The certified optimum is represented by sup - epsilon.
            m_model = res.m_model;
            inf_eps lo(rational(0), inf_rational(res.m_sup_lower, rational(-1)));
            set_best(idx, lo, is_maximize);
            m_lower[idx].set_exact(lo, res.m_sup);
            m_upper[idx].set_exact(inf_eps(rational(0), inf_rational(res.m_sup_upper, rational(-1))), res.m_sup);
            IF_VERBOSE(1, verbose_stream() << "(optsmt nlsat open supremum " << res.m_sup << ")\n");
            return nlsat_outcome::certified;
        }
        inf_eps v(res.m_lower);
        if (v < lower(idx))
            return nlsat_outcome::no_result;
        m_model = res.m_model;
        set_best(idx, v, is_maximize);
        if (r == l_true) {
            m_lower[idx].set_exact(v, res.m_value);
            m_upper[idx].set_exact(inf_eps(res.m_upper), res.m_value);
            IF_VERBOSE(1, verbose_stream() << "(optsmt nlsat optimum " << res.m_value << ")\n");
            return nlsat_outcome::certified;
        }
        // not closed: report the interval [best, sup] where sup is the
        // supremum proven by nlsat when available, else the caller's bound.
        m_upper[idx] = res.m_has_sup ? std::min(inf_eps(res.m_sup_upper), hi) : hi;
        if (upper(idx) < lower(idx))
            m_upper[idx] = m_lower[idx];
        IF_VERBOSE(1, verbose_stream() << "(optsmt nlsat interval [" << lower(idx) << ", " << upper(idx) << "])\n");
        return nlsat_outcome::incomplete;
    }

    /**
       \brief Search the interval [m_lower[idx], hi] for the optimum of a
       real-valued objective when the arithmetic solver cannot push it
       (nonlinear constraints, mod, to_int, ...). Both endpoints are sound:
       the lower bound is backed by m_best_model and hi is a certified upper
       bound, not necessarily a refuted threshold. Alternates two queries
       (GOMT F-Sat / F-Close / F-Split):
         - obj > lo: unsat proves lo optimal; sat raises lo to the new model;
         - obj >= (lo + hi)/2: unsat lowers hi; sat raises lo.
       Returns l_true with m_lower = m_upper when the optimum is proven, and
       l_undef with m_lower < m_upper (the remaining gap) when the budget
       m_bisect_rounds is exhausted or the solver gives up.
    */
    lbool optsmt::bisect(unsigned idx, bool is_maximize, inf_eps hi) {
        arith_util arith(m);
        inf_eps lo = lower(idx);
        inf_eps const eps(rational(0), inf_rational(rational(0), rational(1)));
        SASSERT(lo.is_finite() && hi.is_finite() && lo < hi);
        objective_search search(hi);
        IF_VERBOSE(2, verbose_stream() << "(optsmt bisect [" << lo << ", " << hi << "])\n");
        bool strict = true;
        lbool result = l_undef;
        for (unsigned rounds = 0; m.inc() && lo < hi && rounds < m_bisect_rounds; ++rounds, strict = !strict) {
            inf_eps bound = lo + eps;
            if (!strict) {
                bound = lo + hi;
                bound /= rational(2);
            }
            lbool is_sat;
            {
                solver::scoped_push scope(*m_s);
                search.enter_region(strict ? lo : bound);
                m_s->assert_expr(m_s->mk_ge(idx, bound));
                is_sat = m_s->check_sat(0, nullptr);
                TRACE(opt, tout << "bisect " << (strict ? "strict " : "mid ") << bound << " " << is_sat << "\n";);
                if (is_sat == l_true)
                    m_s->get_model(m_model);
            }
            search.leave_region(is_sat);
            hi = search.upper();
            m_upper[idx] = hi;
            // Publish an improved model only after removing the probe bound.
            if (is_sat == l_true) {
                if (!m_model)
                    break;
                // Prefer the model's own value of the objective; if the model
                // is algebraic (nlsat) and the value is not a numeral, fall
                // back to the bound the model is known to satisfy.
                inf_eps v = bound;
                rational q;
                expr_ref val = (*m_model)(m_objs.get(idx));
                if (arith.is_numeral(val, q) && inf_eps(q) >= bound)
                    v = inf_eps(q);
                else if (strict)
                    continue;  // no exact value to record; the midpoint query makes progress
                if (v > lo) {
                    lo = v;
                    set_best(idx, v, is_maximize);
                }
            }
            else if (is_sat == l_false) {
                if (strict) {
                    // no model is strictly better than the best model.
                    result = l_true;
                    break;
                }
            }
            else
                break;
        }
        m_lower[idx] = lo;
        m_upper[idx] = hi;
        if (result == l_true || lo == hi) 
            return l_true;
        IF_VERBOSE(1, verbose_stream() << "(optsmt interval [" << lo << ", " << hi << "])\n");
        return l_undef;
    }

    bool optsmt::can_increment_delta(vector<inf_eps> const& previous_lower, unsigned i) {
        arith_util arith(m);
        inf_eps max_delta;
        if (lower(i) < upper(i) && arith.is_int(m_objs.get(i))) {
            inf_eps delta = lower(i) - previous_lower[i];
            if (m_lower[i].is_finite() && delta > max_delta) {
                return true;
            }
        }
        return false;
    }

    lbool optsmt::symba_opt() {

        smt::theory_opt& opt = m_s->get_optimizer();

        if (typeid(smt::theory_inf_arith) != typeid(opt)) {
            m_s->set_reason_unknown("symba optimization requires theory_inf_arith");
            return l_undef;
        }


        expr_ref_vector ors(m), disj(m);
        expr_ref fml(m), bound(m.mk_true(), m), tmp(m);
        expr* vars[1];
        {
            for (unsigned i = 0; i < m_upper.size(); ++i) 
                ors.push_back(m_s->mk_ge(i, upper(i)));
                        
            fml = mk_or(ors);
            tmp = m.mk_fresh_const("b", m.mk_bool_sort());
            fml = m.mk_implies(tmp, fml);
            vars[0] = tmp;
            lbool is_sat = l_true;

            solver::scoped_push _push(*m_s);
            while (m.inc()) {
                m_s->assert_expr(fml);
                TRACE(opt, tout << fml << "\n";);
                is_sat = m_s->check_sat(1,vars);
                if (is_sat == l_true) {
                    disj.reset();
                    if (!m_s->maximize_objectives1(disj)) 
                        return l_undef;
                    m_s->get_model(m_model);       
                    m_s->get_labels(m_labels);            
                    for (unsigned i = 0; i < ors.size(); ++i) {
                        if (m_model->is_true(ors.get(i))) {       
                            m_lower[i] = m_upper[i];
                            ors[i]  = m.mk_false();
                            disj[i] = m.mk_false();
                        }
                    }
                    set_max(m_lower, m_s->get_objective_values(), disj);
                    fml = mk_or(ors);
                    tmp = m.mk_fresh_const("b", m.mk_bool_sort());
                    fml = m.mk_implies(tmp, fml);
                    vars[0] = tmp;
                }
                else if (is_sat == l_undef) {
                    return l_undef;
                }
                else {
                    break;
                }
            }
        }
        bound = mk_or(m_lower_fmls);
        m_s->assert_expr(bound);
        
        if (!m.inc()) {
            return l_undef;
        }
        return geometric_opt();
    }

    void optsmt::update_lower_lex(unsigned idx, inf_eps const& v, bool is_maximize) {
        TRACE(opt, tout << v << " lower: " << lower(idx) << "\n";);
        if (v > lower(idx)) {
            m_lower[idx] = v;
            IF_VERBOSE(1, 
                       if (is_maximize) 
                           verbose_stream() << "(optsmt lower bound: " << v << ")\n";
                       else
                           verbose_stream() << "(optsmt upper bound: " << (-v) << ")\n";
                       );            
            for (unsigned i = idx+1; i < m_vars.size(); ++i) {
                m_lower[i] = m_s->saved_objective_value(i);
            }
            TRACE(opt, tout << "update best model " << *m_model << "\n";);
            publish_best_model();
        }
    }

    void optsmt::update_lower(unsigned idx, inf_eps const& v) {
        TRACE(opt, tout << "v" << idx << " >= " << v << "\n";);
        m_lower_fmls[idx] = m_s->mk_ge(idx, v);
        m_lower[idx].update_rational_bound(v);
    }

    void optsmt::update_upper(unsigned idx, inf_eps const& v) {
        TRACE(opt, tout << "v" << idx << " <= " << v << "\n";);
        // A model callback can invalidate this bound by resetting it to
        // infinity; that is not a refinement of its exact finite endpoint.
        if (v.is_finite())
            m_upper[idx].update_rational_bound(v);
        else
            m_upper[idx] = v;
    }

    std::ostream& operator<<(std::ostream& out, vector<objective_value> const& vs) {
        for (unsigned i = 0; i < vs.size(); ++i) {
            out << vs[i].rational_bound() << " ";
        }
        return out;
    }

    expr_ref optsmt::update_lower() {
        expr_ref_vector disj(m);
        m_s->get_model(m_model);
        m_s->get_labels(m_labels);
        if (!m_s->maximize_objectives1(disj))
            return expr_ref(m.mk_true(), m);
        set_max(m_lower, m_s->get_objective_values(), disj);
        TRACE(opt, model_pp(tout << m_lower << "\n", *m_model););
        IF_VERBOSE(2, verbose_stream() << "(optsmt.lower " << m_lower << ")\n";);
        return mk_or(disj);
    }

    void optsmt::setup(opt_solver& solver) {
        m_s = &solver;
        solver.reset_objectives();
        m_vars.reset();

        // force base level
        {
            solver::scoped_push _push(solver);
        }

        for (unsigned i = 0; i < m_objs.size(); ++i) {
            smt::theory_var v = solver.add_objective(m_objs.get(i));
            if (v == smt::null_theory_var) {
                std::ostringstream out;
                out << "Objective function '" << mk_pp(m_objs.get(i), m) << "' is not supported";
                throw default_exception(out.str());
            }
            m_vars.push_back(v);
        }            
    }

    lbool optsmt::lex(unsigned obj_index, bool is_maximize) {
        TRACE(opt, tout << "optsmt:lex\n";);
        m_context.get_base_model(m_best_model);
        solver::scoped_push _push(*m_s);
        SASSERT(obj_index < m_vars.size());
        if (is_maximize && m_optsmt_engine == symbol("symba"))
            return symba_opt();

        // A fresh search must not report a previous exact value for this objective.
        m_lower[obj_index].reset_exact();
        m_upper[obj_index].reset_exact();
        // Add constraints so optimizing this objective cannot worsen earlier
        // objectives' results. Other variable values remain free to change.
        for (unsigned i = 0; i < obj_index; ++i)
            commit_assignment(i);

        lbool result = geometric_search(obj_index, is_maximize);
        // Later objectives must be reconsidered only after this one succeeds.
        if (result == l_true)
            for (unsigned i = obj_index + 1; i < m_lower.size(); ++i)
                m_lower[i] = inf_eps(rational(-1), inf_rational(0));
        return result;
    }

    // Under the current assertions, does objective idx still have the
    // supremum +oo? All objective bounds are restored afterwards. On a
    // positive answer, witness is a feasible model from the search.
    lbool optsmt::reaches_infinity(unsigned idx, model_ref& witness) {
        vector<objective_value> lower(m_lower), upper(m_upper);
        model_ref best = m_best_model, model = m_model;
        m_lower[idx] = inf_eps(rational(-1), inf_rational(0));
        m_upper[idx] = inf_eps(rational(1), inf_rational(0));
        // geometric_search reports infeasibility only while m_model is null.
        m_model = nullptr;
        lbool r;
        {
            solver::scoped_push _push(*m_s);
            // The objective is known to be unbounded over F; one climb round
            // precedes the proof attempt.
            r = geometric_search(idx, true, 1);
        }
        bool unbounded = r == l_true && !m_upper[idx].is_finite();
        if (unbounded)
            witness = m_best_model;
        m_lower = lower;
        m_upper = upper;
        m_best_model = best;
        m_model = model;
        if (r == l_undef)
            return l_undef;
        return unbounded ? l_true : l_false;
    }

    bool optsmt::model_objective_value(unsigned idx, model_ref& mdl, rational& value) {
        if (!mdl)
            return false;
        arith_util arith(m);
        expr_ref val = (*mdl)(m_objs.get(idx));
        return model_value_bound(arith, val, true, value);
    }

    /**
       Objective u = unbounded_index has the supremum +oo over the current
       assertions F. Objective f = obj_index then takes the value

           L = lim_{c -> +oo} sup { f : F and u > c }.

       For a rational t: if F and f >= t keeps u unbounded, then t <= L;
       if F and f > t leaves u bounded or infeasible, then L <= t, because
       beyond some value of u no feasible point has f > t. A threshold with
       both facts is L itself, attained at every feasible point with f = t
       and arbitrarily large u. Thresholds come from model values and from
       the arithmetic solver's hints, and every threshold is checked by the
       unbounded objective's own search. A limit of +oo or -oo is not
       certified: a streak of trials in the same direction ends the search.
       Without a deciding threshold the objective keeps the sound interval
       [lo, hi] and the result stays sat.
    */
    lbool optsmt::lex_after_unbounded(unsigned unbounded_index, unsigned idx, bool is_maximize) {
        TRACE(opt, tout << "optsmt:lex after unbounded " << unbounded_index << " objective " << idx << "\n";);
        if (!m_lex_after_unbounded)
            return l_true;
        arith_util arith(m);
        app* obj = m_objs.get(idx);
        bool is_int = arith.is_int(obj);
        inf_eps const minus_infty(rational(-1), inf_rational(0));
        inf_eps const infty(rational(1), inf_rational(0));
        inf_eps lo = minus_infty, hi = infty;
        std::optional<rational> candidate;
        // f > lo keeps u unbounded, so lo is not the limit unless hi closes in.
        bool above_lo_feasible = false;
        model_ref witness, best;
        geometric_step step;
        scoped_pushes scopes(*m_s);
        lbool result = l_undef;
        bool decided = false;
        unsigned const max_one_sided_rounds = 16;
        unsigned one_sided_rounds = 0;
        for (unsigned rounds = 0; m.inc() && rounds < m_bisect_rounds; ++rounds) {
            if (is_int && lo.is_finite() && hi.is_finite() && hi.get_rational() <= lo.get_rational() + 1) {
                decided = true;
                break;
            }
            if (best && (!lo.is_finite() || !hi.is_finite()) && ++one_sided_rounds > max_one_sided_rounds)
                break;
            bool has_trial = true, strict = false;
            rational t;
            if (candidate && inf_eps(*candidate) > lo && inf_eps(*candidate) < hi)
                t = *candidate;
            else if (!best)
                has_trial = false;
            else if (!lo.is_finite())
                t = hi.get_rational() - step.value();
            else if (!hi.is_finite())
                t = lo.get_rational() + step.value();
            else if (!is_int && !above_lo_feasible) {
                t = lo.get_rational();
                strict = true;
            }
            else {
                t = (lo.get_rational() + hi.get_rational()) / rational(2);
                if (is_int)
                    t = floor(t);
            }
            if (is_int && has_trial)
                t = ceil(t);
            candidate.reset();
            scopes.push();
            if (has_trial) {
                expr_ref num(arith.mk_numeral(t, is_int), m);
                expr_ref bound(strict ? arith.mk_gt(obj, num) : arith.mk_ge(obj, num), m);
                TRACE(opt, tout << "trial " << bound << "\n";);
                m_s->assert_expr(bound);
            }
            lbool r = reaches_infinity(unbounded_index, witness);
            inf_eps hint;
            bool has_hint = false;
            if (r == l_true && m_s->check_sat(0, nullptr) == l_true) {
                // The largest value of f the arithmetic solver sees under the
                // trial is only a threshold to check, not a certified bound.
                expr_ref blocker(m);
                auto hint_result = m_s->maximize_objective(idx, blocker, false, false);
                if (hint_result.bound_valid && hint_result.hint.is_finite()) {
                    hint = hint_result.hint;
                    has_hint = true;
                }
            }
            scopes.pop();
            if (r == l_undef) {
                result = l_undef;
                break;
            }
            result = l_true;
            if (r == l_true) {
                best = witness;
                if (has_trial && strict)
                    above_lo_feasible = true;
                else if (has_trial) {
                    lo = inf_eps(t);
                    above_lo_feasible = false;
                }
                rational v;
                if (model_objective_value(idx, best, v) && inf_eps(v) > lo)
                    candidate = v;
                if (has_hint && hint > lo && (!candidate || hint.get_rational() > *candidate))
                    candidate = hint.get_rational();
                step.update(true);
                IF_VERBOSE(2, verbose_stream() << "(optsmt.lex-after-unbounded :lower " << lo << " :upper " << hi << ")\n";);
            }
            else {
                if (!has_trial)
                    break;
                if (strict) {
                    decided = true;
                    break;
                }
                hi = inf_eps(t);
                step.reset();
                IF_VERBOSE(2, verbose_stream() << "(optsmt.lex-after-unbounded :lower " << lo << " :upper " << hi << ")\n";);
            }
        }
        scopes.reset();
        m_lower[idx] = lo;
        m_upper[idx] = decided ? lo : hi;
        IF_VERBOSE(1, verbose_stream() << "(optsmt " << (is_maximize ? "lower" : "upper") << " bound after unbounded objective: "
                   << (is_maximize ? lo : -lo) << ")\n";);
        if (best) {
            m_model = best;
            publish_best_model();
        }
        return result;
    }

    /**
       Takes solver with hard constraints added.
       Returns an optimal assignment to objective functions.
    */
    lbool optsmt::box() {
        lbool is_sat = l_true;        
        if (m_vars.empty()) {
            return is_sat;
        }
        // In box mode, optimize each objective independently.
        // Each objective gets its own push/pop scope so that bounds
        // from one objective do not constrain another.
        // Note: geometric_search is used unconditionally here, even when
        // m_optsmt_engine is "symba", because symba_opt and geometric_opt
        // optimize all objectives jointly, violating box mode semantics.
        //
        m_context.get_base_model(m_best_model);
        for (unsigned i = 0; i < m_vars.size() && m.inc(); ++i) {
            // Reset bounds for objective i so that update_lower_lex
            // contamination from earlier objectives does not affect it.
            m_lower[i] = inf_eps(rational(-1), inf_rational(0));
            m_upper[i] = inf_eps(rational(1), inf_rational(0));
            solver::scoped_push _push(*m_s);
            is_sat = geometric_search(i, true);
            if (is_sat == l_undef)
                return l_undef;
            if (is_sat == l_false)
                return l_false;
            m_models.set(i, m_best_model.get());
        }
        return l_true;
    }


    objective_value optsmt::get_lower(unsigned i) const {
        if (i >= m_lower.size())
            return objective_value(m);
        return m_lower[i];
    }

    objective_value optsmt::get_upper(unsigned i) const {
        if (i >= m_upper.size())
            return objective_value(m);
        return m_upper[i];
    }

    void optsmt::get_model(model_ref& mdl, svector<symbol> & labels) {        
        mdl = m_best_model.get();
        TRACE(opt, tout << *mdl << "\n";);
        labels = m_labels;
    }

    // force lower_bound(i) <= objective_value(i)    
    void optsmt::commit_assignment(unsigned i) {
        inf_eps lo = lower(i);
        TRACE(opt, tout << "set lower bound of " << mk_pp(m_objs.get(i), m) << " to: " << lo << "\n";
              tout << lower(i) << ":" << upper(i) << "\n";);
        expr* e = m_lower[i].exact_finite();
        SASSERT(!has_open_bound(i));
        arith_util arith(m);
        if (e && arith.is_irrational_algebraic_numeral(e)) {
            // Keep the defining equation and its isolating interval in
            // rational arithmetic; the SMT core cannot use a root numeral directly.
            expr_ref constraint = mk_algebraic_eq(m, m_objs.get(i), e);
            m_s->assert_expr(constraint);
            TRACE(opt, tout << "commit exact " << mk_pp(e, m) << ": " << constraint << "\n";);
        }
        else if (lo.is_finite()) {
            // Only assert bounds for bounded objectives
            m_s->assert_expr(m_s->mk_ge(i, lo));
        }
    }

    unsigned optsmt::add(app* t) {
        expr_ref t1(t, m), t2(m);
        th_rewriter rw(m);
        rw(t1, t2);
        SASSERT(is_app(t2));
        m_objs.push_back(to_app(t2));
        m_lower.push_back(objective_value(m, inf_eps(rational(-1),inf_rational(0))));
        m_upper.push_back(objective_value(m, inf_eps(rational(1), inf_rational(0))));
        m_lower_fmls.push_back(m.mk_true());
        m_models.push_back(nullptr);
        return m_objs.size()-1;
    }

    void optsmt::updt_params(params_ref& p) {
        opt_params _p(p);
        m_optsmt_engine = _p.optsmt_engine();        
        m_bisect_rounds = _p.optsmt_bisect_rounds();
        m_optsmt_nlsat = _p.optsmt_nlsat();
        m_dual_bounds = _p.optsmt_dual_bounds();
        m_nlsat_supremum_rlimit = _p.optsmt_nlsat_supremum_rlimit();
        m_lex_after_unbounded = _p.optsmt_lex_after_unbounded();
    }

    void optsmt::reset() {
        m_lower.reset();
        m_upper.reset();
        m_objs.reset();
        m_vars.reset();
        m_model.reset();       
        m_best_model = nullptr; 
        m_models.reset();
        m_lower_fmls.reset();
        m_s = nullptr;
    }
}