Program Listing for File cutting_plane.hpp

↰ Return to documentation for file (ellcpp/cutting_plane.hpp)

// -*- coding: utf-8 -*-
#pragma once

#include "cut_config.hpp"
#include "half_nonnegative.hpp"
#include <cassert>
#include <cmath>
#include <tuple>


template <typename Oracle, typename Space>
auto cutting_plane_feas(
    Oracle&& Omega, Space&& S, const Options& options = Options()) -> CInfo
{
    auto feasible = false;
    auto status = CUTStatus::success;

    auto niter = 0U;
    while (++niter != options.max_it)
    {
        const auto cut = Omega(S.xc()); // query the oracle at S.xc()
        if (!cut)
        { // feasible sol'n obtained
            feasible = true;
            break;
        }
        const auto [cutstatus, tsq] = S.update(*cut); // update S
        if (cutstatus != CUTStatus::success)
        {
            status = cutstatus;
            break;
        }
        if (tsq < options.tol)
        { // no more
            status = CUTStatus::smallenough;
            break;
        }
    }
    return {feasible, niter, status};
}

template <typename Oracle, typename Space, typename opt_type>
auto cutting_plane_dc(
    Oracle&& Omega, Space&& S, opt_type&& t, const Options& options = Options())
{
    const auto t_orig = t;
    decltype(S.xc()) x_best;
    auto status = CUTStatus::success;

    auto niter = 0U;
    while (++niter != options.max_it)
    {
        const auto [cut, shrunk] = Omega(S.xc(), t);
        if (shrunk)
        { // best t obtained
            x_best = S.xc();
        }
        const auto [cutstatus, tsq] = S.update(cut);
        if (cutstatus != CUTStatus::success) // ???
        {
            status = cutstatus;
            break;
        }
        if (tsq < options.tol)
        { // no more
            status = CUTStatus::smallenough;
            break;
        }
    }
    return std::make_tuple(
        std::move(x_best), CInfo {t != t_orig, niter, status});
} // END

// #include <boost/numeric/ublas/symmetric.hpp>
// namespace bnu = boost::numeric::ublas;
// #include <xtensor-blas/xlinalg.hpp>
// #include <xtensor/xarray.hpp>

template <typename Oracle, typename Space, typename opt_type>
auto cutting_plane_q(
    Oracle&& Omega, Space&& S, opt_type&& t, const Options& options = Options())
{
    const auto t_orig = t;
    decltype(S.xc()) x_best;
    auto status = CUTStatus::nosoln; // note!!!

    auto niter = 0U;
    while (++niter != options.max_it)
    {
        auto retry = (status == CUTStatus::noeffect);
        const auto [cut, x0, shrunk, more_alt] = Omega(S.xc(), t, retry);
        if (shrunk)
        { // best t obtained
            // t = t1;
            x_best = x0;
        }
        const auto [cutstatus, tsq] = S.update(cut);
        if (cutstatus == CUTStatus::noeffect)
        {
            if (!more_alt)
            {
                break; // no more alternative cut
            }
        }
        if (cutstatus == CUTStatus::nosoln)
        {
            status = cutstatus;
            break;
        }
        if (tsq < options.tol)
        {
            status = CUTStatus::smallenough;
            break;
        }
    }
    return std::make_tuple(
        std::move(x_best), CInfo {t != t_orig, niter, status});
} // END

template <typename Oracle, typename Space>
auto bsearch(Oracle&& Omega, Space&& I, const Options& options = Options())
    -> CInfo
{
    // assume monotone
    // auto& [lower, upper] = I;
    auto& lower = I.first;
    auto& upper = I.second;
    assert(lower <= upper);
    const auto u_orig = upper;
    auto niter = 0U;
    auto status = CUTStatus::success;

    for (; niter != options.max_it; ++niter)
    {
        auto tau = algo::half_nonnegative(upper - lower);
        if (tau < options.tol)
        {
            status = CUTStatus::smallenough;
            break;
        }

        auto t = lower; // l may be `int` or `Fraction`
        t += tau;
        if (Omega(t))
        { // feasible sol'n obtained
            upper = t;
        }
        else
        {
            lower = t;
        }
    }
    return {upper != u_orig, niter + 1, status};
}

template <typename Oracle, typename Space> //
class bsearch_adaptor
{
  private:
    Oracle& _P;
    Space& _S;
    const Options _options;

  public:
    bsearch_adaptor(Oracle& P, Space& S)
        : bsearch_adaptor {P, S, Options()}
    {
    }

    bsearch_adaptor(Oracle& P, Space& S, const Options& options)
        : _P {P}
        , _S {S}
        , _options {options}
    {
    }

    [[nodiscard]] auto x_best() const
    {
        return this->_S.xc();
    }

    template <typename opt_type>
    auto operator()(const opt_type& t) -> bool
    {
        Space S = this->_S.copy();
        this->_P.update(t);
        const auto ell_info = cutting_plane_feas(this->_P, S, this->_options);
        if (ell_info.feasible)
        {
            this->_S.set_xc(S.xc());
        }
        return ell_info.feasible;
    }
};