cutting_plane.hpp

Go to the documentation of this file.

 
#pragma once
 
#include <cassert>
#include <cmath>
#include <tuple>
#include <type_traits>
 
#include "ell_config.hpp"
#include "half_nonnegative.hpp"
 
template <typename S> using CuttingPlaneArrayType = typename std::remove_reference_t<S>::ArrayType;
 
template <typename T> inline auto invalid_value() -> T
    requires std::is_floating_point_v<T>
{
    return std::nan("1");
}
 
template <typename T> inline auto invalid_value() -> T
    requires(!std::is_floating_point_v<T>)
{
    return T{};
}
 
template <typename O, typename S>
    requires OracleFeas<O, typename S::ArrayType> && SearchSpace<S>
inline auto cutting_plane_feas(O& omega, S& space, const Options& options = Options())
    -> std::tuple<CuttingPlaneArrayType<S>, size_t> {
    for (auto niter = 0U; niter != options.max_iters; ++niter) {
        const auto cut = omega.assess_feas(space.xc());
        if (!cut) {  // feasible sol'n obtained
            return {space.xc(), niter};
        }
        const auto status = space.update_bias_cut(*cut);  // update space
        if (status != CutStatus::Success || space.tsq() < options.tolerance) {
            auto res = invalid_value<CuttingPlaneArrayType<S>>();
            return {std::move(res), niter};
        }
    }
    auto res = invalid_value<CuttingPlaneArrayType<S>>();
    return {std::move(res), options.max_iters};
}
 
template <typename O, typename S, typename N>
    requires OracleOptim<O, typename S::ArrayType, N> && SearchSpace<S>
inline auto cutting_plane_optim(O& omega, S& space, N& gamma, const Options& options = Options())
    -> std::tuple<CuttingPlaneArrayType<S>, size_t> {
    auto x_best = invalid_value<CuttingPlaneArrayType<S>>();
    for (auto niter = 0U; niter < options.max_iters; ++niter) {
        const auto _result1 = omega.assess_optim(space.xc(), gamma);
        const auto& cut = std::get<0>(_result1);
        const auto& shrunk = std::get<1>(_result1);
        const auto status = [&]() {
            if (shrunk) {  // best gamma obtained
                x_best = space.xc();
                return space.update_central_cut(cut);  // should update_central_cut
            }
            return space.update_bias_cut(cut);
        }();
        if (status != CutStatus::Success || space.tsq() < options.tolerance) {  // no more
            return {std::move(x_best), niter};
        }
    }
    return {std::move(x_best), options.max_iters};
}  // END
 
template <typename A> class OptimQState {
  public:
    enum class Result { Continue, NoSoln, NoMoreAlt };
 
  private:
    A _x_best;
    bool _retry = false;
 
  public:
    explicit OptimQState(A invalid) : _x_best{std::move(invalid)} {}
 
    auto x_best() const -> const A& { return this->_x_best; }
 
    auto x_best() -> A& { return this->_x_best; }
 
    auto retry() const -> bool { return this->_retry; }
 
    void on_shrunk(A x) {
        this->_x_best = std::move(x);
        this->_retry = false;
    }
 
    auto on_update(const CutStatus status, const bool more_alt) -> Result {
        switch (status) {
            case CutStatus::Success:
                this->_retry = false;
                return Result::Continue;
            case CutStatus::NoSoln:
                return Result::NoSoln;
            case CutStatus::NoEffect:
                if (more_alt) {
                    this->_retry = true;
                    return Result::Continue;
                }
                return Result::NoMoreAlt;
            case CutStatus::Unknown:
                return Result::Continue;
        }
        return Result::Continue;
    }
};
 
template <typename O, typename S, typename N>
    requires OracleOptimQ<O, typename S::ArrayType, N> && SearchSpace<S>
inline auto cutting_plane_optim_q(O& omega, S& space_q, N& gamma,
                                  const Options& options = Options())
    -> std::tuple<CuttingPlaneArrayType<S>, size_t> {
    using A = CuttingPlaneArrayType<S>;
    OptimQState<A> state{invalid_value<A>()};
 
    for (auto niter = 0U; niter < options.max_iters; ++niter) {
        const auto result1 = omega.assess_optim_q(space_q.xc(), gamma, state.retry());
        const auto& cut = std::get<0>(result1);
        const auto& shrunk = std::get<1>(result1);
        if (shrunk) {  // best gamma obtained
            state.on_shrunk(std::move(std::get<2>(result1)));
        }
        const auto outcome = state.on_update(space_q.update_q(cut), std::get<3>(result1));
        if (outcome == OptimQState<A>::Result::NoSoln) {
            return {std::move(state.x_best()), niter};
        }
        if (outcome == OptimQState<A>::Result::NoMoreAlt) {
            break;  // no more alternative cut
        }
        if (space_q.tsq() < options.tolerance) {  // no more
            return {std::move(state.x_best()), niter};
        }
    }
    return {std::move(state.x_best()), options.max_iters};
}  // END
 
template <typename Oracle, typename Space>  //
class BSearchAdaptor {
    using ArrayType = typename Space::ArrayType;
 
    Oracle* _omega;
    Space* _space;
    Options _options;
 
  public:
    BSearchAdaptor(Oracle& omega, Space& space) : BSearchAdaptor{omega, space, Options()} {}
 
    BSearchAdaptor(Oracle& omega, Space& space, const Options& options)
        : _omega{&omega}, _space{&space}, _options{options} {}
 
    auto x_best() const -> ArrayType { return this->_space->xc(); }
 
    template <typename Num> auto assess_bs(Num& gamma) -> bool {
        Space space = this->_space->copy();  // copy
        this->_omega->update(gamma);
        const auto result = cutting_plane_feas(*this->_omega, space, this->_options);
        auto x_feas = std::get<0>(result);
        if (x_feas.size() != 0U) {
            this->_space->set_xc(x_feas);
            return true;
        }
        return false;
    }
};
 
template <typename O, typename T>
    requires OracleBS<O, T>
inline auto bsearch(O& omega, const std::pair<T, T>& intvl, const Options& options = Options())
    -> std::tuple<T, size_t> {
    // assume monotone
    // auto& [lower, upper] = intvl;
    auto lower = intvl.first;
    auto upper = intvl.second;
    assert(lower <= upper);
 
    for (auto niter = 0U; niter < options.max_iters; ++niter) {
        auto tau = algo::half_nonnegative(upper - lower);
        if (tau < options.tolerance) {  // no more
            return {upper, niter};
        }
        auto gamma = lower;  // l may be `int` or `Fraction`
        gamma += tau;
        if (omega.assess_bs(gamma)) {  // feasible sol'n obtained
            upper = gamma;
        } else {
            lower = gamma;
        }
    }
    return {upper, options.max_iters};
}