:github_url: https://github.com/svenevs/exhale-companion


.. _program_listing_file_ellcpp_cutting_plane.hpp:

Program Listing for File cutting_plane.hpp
==========================================

|exhale_lsh| :ref:`Return to documentation for file <file_ellcpp_cutting_plane.hpp>` (``ellcpp/cutting_plane.hpp``)

.. |exhale_lsh| unicode:: U+021B0 .. UPWARDS ARROW WITH TIP LEFTWARDS

.. code-block:: cpp

   // -*- 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;
       }
   };