itm-common.hpp

/* Copyright (C) 2016-2019 INRA
 *
 * Permission is hereby granted, free of charge, to any person obtaining a
 * copy of this software and associated documentation files (the
 * "Software"), to deal in the Software without restriction, including
 * without limitation the rights to use, copy, modify, merge, publish,
 * distribute, sublicense, and/or sell copies of the Software, and to
 * permit persons to whom the Software is furnished to do so, subject to
 * the following conditions:
 *
 * The above copyright notice and this permission notice shall be included
 * in all copies or substantial portions of the Software.
 *
 * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
 * OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF
 * MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND
 * NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE
 * LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION
 * OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION
 * WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.
 */

#ifndef ORG_VLEPROJECT_BARYONYX_SOLVER_ITM_COMMON_HPP
#define ORG_VLEPROJECT_BARYONYX_SOLVER_ITM_COMMON_HPP

#include <baryonyx/core>

#include "bit-array.hpp"
#include "debug.hpp"
#include "private.hpp"
#include "problem.hpp"
#include "sparse-matrix.hpp"
#include "utils.hpp"

#include <algorithm>
#include <chrono>
#include <fstream>
#include <iterator>
#include <numeric>
#include <random>
#include <thread>
#include <tuple>
#include <vector>

#include <fmt/ostream.h>

namespace baryonyx {
namespace itm {

using random_engine = std::default_random_engine;

struct maximize_tag
{};

struct minimize_tag
{};

struct merged_constraint
{
    merged_constraint(std::vector<function_element> elements_,
                      int min_,
                      int max_,
                      int id_)
      : elements(std::move(elements_))
      , min(min_)
      , max(max_)
      , id(id_)
    {}

    std::vector<function_element> elements;
    int min;
    int max;
    int id;
};

std::vector<merged_constraint>
make_merged_constraints(const context_ptr& ctx, const problem& pb);

template<typename Solver, typename Xtype>
bool
is_valid_constraint(const Solver& slv, int k, const Xtype& x)
{
    typename sparse_matrix<int>::const_row_iterator it, et;

    std::tie(it, et) = slv.ap.row(k);
    int v = 0;

    for (; it != et; ++it)
        v += slv.factor(it->value) * x[it->column];

    return slv.bound_min(k) <= v && v <= slv.bound_max(k);
}

template<typename Solver, typename Xtype>
bool
is_valid_solution(const Solver& slv, const Xtype& x)
{
    for (int k = 0; k != slv.m; ++k)
        if (!is_valid_constraint(slv, k, x))
            return false;

    return true;
}

template<typename Solver, typename Xtype>
int
compute_violated_constraints(const Solver& slv,
                             const Xtype& x,
                             std::vector<int>& out)
{
    out.clear();

    for (int k = 0; k != slv.m; ++k)
        if (!is_valid_constraint(slv, k, x))
            out.emplace_back(k);

    return length(out);
}

template<typename iteratorT>
inline void
random_shuffle_unique(iteratorT begin,
                      iteratorT end,
                      random_engine& rng) noexcept
{
    auto ret = begin++;
    for (; begin != end; ++begin) {
        if (ret->value != begin->value) {
            std::shuffle(ret, begin, rng);
            ret = begin;
        }
    }

    std::shuffle(ret, begin, rng);
}

template<typename Mode, typename Iterator>
inline void
calculator_sort(Iterator begin, Iterator end, random_engine& rng)
{
    if (std::distance(begin, end) > 1) {
        std::sort(begin, end, [](const auto& lhs, const auto& rhs) {
            if constexpr (std::is_same_v<Mode, minimize_tag>)
                return lhs.value < rhs.value;
            else
                return rhs.value < lhs.value;
        });

        random_shuffle_unique(begin, end, rng);
    }
}

template<typename Mode, typename Float>
inline bool
stop_iterating(Float value, random_engine& rng) noexcept
{
    if (value == 0) {
        std::bernoulli_distribution d(0.5);
        return d(rng);
    }

    if constexpr (std::is_same_v<Mode, minimize_tag>)
        return value > 0;
    else
        return value < 0;
}

template<typename Mode, typename Float>
constexpr Float
bad_value() noexcept
{
    if constexpr (std::is_same_v<Mode, minimize_tag>)
        return std::numeric_limits<Float>::infinity();
    else
        return -std::numeric_limits<Float>::infinity();
}

template<typename Mode, typename Float>
constexpr bool
is_bad_value(Float value) noexcept
{
    return value == bad_value<Mode, Float>();
}

template<typename Mode, typename Float>
inline bool
stop_iterating(Float value) noexcept
{
    if constexpr (std::is_same_v<Mode, minimize_tag>)
        return value > 0;
    else
        return value < 0;
}

template<typename Mode, typename Float>
inline bool
is_better_solution(Float lhs, Float rhs) noexcept
{
    if constexpr (std::is_same_v<Mode, minimize_tag>)
        return lhs < rhs;
    else
        return lhs > rhs;
}

template<typename Mode, typename Float>
inline bool
init_x(Float cost, int value_if_cost_0) noexcept
{
    if constexpr (std::is_same_v<Mode, minimize_tag>) {
        if (cost < 0)
            return true;

        if (cost == 0)
            return value_if_cost_0;
    } else {
        if (cost > 0)
            return true;

        if (cost == 0)
            return value_if_cost_0;
    }

    return false;
}

template<typename TimePoint>
double
compute_duration(const TimePoint& first, const TimePoint& last) noexcept
{
    namespace sc = std::chrono;

    auto diff = last - first;
    auto dc = sc::duration_cast<sc::duration<double, std::ratio<1>>>(diff);

    return dc.count();
}

inline std::size_t
compute_reduced_costs_vector_size(
  const std::vector<merged_constraint>& csts) noexcept
{
    std::size_t rsizemax = csts[0].elements.size();

    for (std::size_t i = 1, e = csts.size(); i != e; ++i)
        rsizemax = std::max(rsizemax, csts[i].elements.size());

    return rsizemax;
}

/// Test if the bit sign between two reals are differents.
template<typename Float>
inline bool
is_signbit_change(Float lhs, Float rhs) noexcept
{
    return std::signbit(lhs) != std::signbit(rhs);
}

/// The bkmin and bkmax constraint bounds are not equal and can be assigned to
/// -infinity or +infinity. We have to scan the r vector and search a value j
/// such as b(0, k) <= Sum A(k, R[j]) < b(1, k).
///
/// @return True if the sign of the pi vector change during the affect
///     operation.
template<typename Solver, typename Xtype, typename Iterator, typename Float>
bool
affect(Solver& slv,
       Xtype& x,
       Iterator it,
       int k,
       int selected,
       int r_size,
       const Float kappa,
       const Float delta)
{
    constexpr Float one{ 1 };
    constexpr Float two{ 2 };
    constexpr Float middle{ (two + one) / two };

    const auto old_pi = slv.pi[k];

    auto d = delta;

    if (selected < 0) {
        // slv.pi[k] += slv.R[0].value / two;
        d += (kappa / (one - kappa)) * (slv.R[0].value / two);

        for (int i = 0; i != r_size; ++i) {
            auto var = it + slv.R[i].id;

            if (slv.R[i].is_negative_factor()) {
                x.set(var->column);
                slv.P[var->value] += d;
            } else {
                x.unset(var->column);
                slv.P[var->value] -= d;
            }
        }
    } else if (selected + 1 >= r_size) {
        // slv.pi[k] += slv.R[selected].value * middle;
        d += (kappa / (one - kappa)) * (slv.R[selected].value * middle);

        for (int i = 0; i != r_size; ++i) {
            auto var = it + slv.R[i].id;

            if (slv.R[i].is_negative_factor()) {
                x.unset(var->column);
                slv.P[var->value] -= d;
            } else {
                x.set(var->column);
                slv.P[var->value] += d;
            }
        }
    } else {
        slv.pi[k] +=
          ((slv.R[selected].value + slv.R[selected + 1].value) / two);
        d += (kappa / (one - kappa)) *
             (slv.R[selected + 1].value - slv.R[selected].value);

        int i = 0;
        for (; i <= selected; ++i) {
            auto var = it + slv.R[i].id;

            if (slv.R[i].is_negative_factor()) {
                x.unset(var->column);
                slv.P[var->value] -= d;
            } else {
                x.set(var->column);
                slv.P[var->value] += d;
            }
        }

        for (; i != r_size; ++i) {
            auto var = it + slv.R[i].id;

            if (slv.R[i].is_negative_factor()) {
                x.set(var->column);
                slv.P[var->value] += d;
            } else {
                x.unset(var->column);
                slv.P[var->value] -= d;
            }
        }
    }

    // TODO job: develops is_valid_constraint for all the solvers
    bx_expects(is_valid_constraint(slv, k, x));

    return is_signbit_change(old_pi, slv.pi[k]);
}

namespace detail {

inline int
constraint(std::vector<int>::const_iterator it)
{
    return *it;
}

inline int
constraint(std::vector<int>::const_reverse_iterator it)
{
    return *it;
}

inline int
constraint(std::vector<std::pair<int, int>>::const_iterator it)
{
    return it->first;
}
}

template<typename Iterator>
int
constraint(Iterator it)
{
    return detail::constraint(it);
}

template<typename Float, typename Mode>
struct best_solution_recorder;

template<typename Solver, typename Float, typename Mode>
class solver_initializer
{
    std::bernoulli_distribution dist;
    std::bernoulli_distribution toss_up;
    bool use_cycle{ true };

    enum class single_automaton : std::uint8_t
    {
        init,
        improve_x_1,
        improve_x_2,
        improve_x_3
    };

    enum class cycle_automaton : std::uint8_t
    {
        bastert,
        bastert_crossover,
        pessimistic_solve,
        pessimistic_crossover,
        optimistic_solve,
        optimistic_crossover
    };

    std::array<unsigned, 4> cycle_counter = { 0, 0, 0, 0 };

    enum class init_automaton : std::uint8_t
    {
        random,
        current_best,
        best_recorder
    };

    static inline const std::string_view single_automaton_string[] =
      { "init", "improve_x_1", "improve_x_2", "improve_x_3" };

    static inline const std::string_view cycle_automaton_string[] = {
        "bastert",           "bastert_crossover",
        "pessimistic_solve", "pessimistic_solve_crossver",
        "optimistic_solve",  "optimistic_solve_crossover",
    };

    static inline const std::string_view init_automaton_string[] =
      { "random", "current_best", "best_recorder" };

    single_automaton single{ single_automaton::init };
    cycle_automaton cycle{ cycle_automaton::pessimistic_solve };
    init_automaton base{ init_automaton::random };

    void init_crossover(
      Solver& slv,
      bit_array& x,
      const best_solution_recorder<Float, Mode>& best_recorder)
    {
        to_log(stdout, 6u, "initializer: crossover {}\n", cycle_counter[3]);
        ++cycle_counter[3];

        const int min = 1;
        const int max = slv.n >= 2 ? slv.n - 1 : 1;

        std::uniform_int_distribution<int> uniform_dist(min, max);
        const int split_at = uniform_dist(slv.rng);

        // TODO Need to improve this part to avoid bit-per-bit copy using the
        // std::uint64_t underlying type.

        if (toss_up(slv.rng)) {
            if (toss_up(slv.rng)) {
                best_recorder.copy_any_solution(x, slv.rng, 0, split_at);
                best_recorder.copy_best_solution(x, split_at, slv.n);
            } else {
                best_recorder.copy_best_solution(x, 0, split_at);
                best_recorder.copy_any_solution(x, slv.rng, split_at, slv.n);
            }
        } else {
            best_recorder.copy_any_solution(x, slv.rng, 0, split_at);
            best_recorder.copy_any_solution(x, slv.rng, split_at, slv.n);
        }
    }

    void init_bastert(Solver& slv, bit_array& x) noexcept
    {
        to_log(stdout, 6u, "- initializer: bastert {}\n", cycle_counter[0]);
        ++cycle_counter[0];

        const int value_if_cost_0 = 1;

        for (int i = 0; i != slv.n; ++i) {
            if (dist(slv.rng)) {
                if (init_x<Mode>(slv.c(i, x), value_if_cost_0))
                    x.set(i);
                else
                    x.unset(i);
            }
        }
    }

    void init_pessimistic_solve(Solver& slv, bit_array& x) noexcept
    {
        to_log(
          stdout, 6u, "- initializer: pessimistic {}\n", cycle_counter[1]);
        ++cycle_counter[1];

        for (int k = 0; k != slv.m; ++k) {
            if (!dist(slv.rng))
                continue;

            auto [begin, end] = slv.ap.row(k);
            int r_size = 0;

            for (auto it = begin; it != end; ++it, ++r_size) {
                slv.R[r_size].value = slv.c(it->column, x);
                slv.R[r_size].id = r_size;
            }

            std::shuffle(slv.R.get(), slv.R.get() + r_size, slv.rng);

            std::sort(slv.R.get(),
                      slv.R.get() + r_size,
                      [](const auto& lhs, const auto& rhs) {
                          if constexpr (std::is_same_v<Mode, minimize_tag>)
                              return lhs.value < rhs.value;
                          else
                              return rhs.value < lhs.value;
                      });

            int sum = 0;
            int best = -2;

            for (int i = -1; i < r_size; ++i) {
                if (slv.bound_min(k) <= sum && sum <= slv.bound_max(k)) {
                    best = i;
                    break;
                }

                auto var = begin + slv.R[i].id;
                sum += slv.factor(var->value);
            }

            int i = 0;
            for (; i <= best; ++i) {
                auto var = begin + slv.R[i].id;
                x.set(var->column);
            }

            for (; i != r_size; ++i) {
                auto var = begin + slv.R[i].id;
                x.unset(var->column);
            }
        }
    }

    void init_optimistic_solve(Solver& slv, bit_array& x) noexcept
    {
        to_log(stdout, 6u, "- initializer: optimistic {}\n", cycle_counter[2]);
        ++cycle_counter[2];

        for (int k = 0; k != slv.m; ++k) {
            if (!dist(slv.rng))
                continue;

            auto [begin, end] = slv.ap.row(k);
            int r_size = 0;

            for (auto it = begin; it != end; ++it, ++r_size) {
                slv.R[r_size].value = slv.c(it->column, x);
                slv.R[r_size].id = r_size;
            }

            std::shuffle(slv.R.get(), slv.R.get() + r_size, slv.rng);

            std::sort(slv.R.get(),
                      slv.R.get() + r_size,
                      [](const auto& lhs, const auto& rhs) {
                          if constexpr (std::is_same_v<Mode, minimize_tag>)
                              return lhs.value < rhs.value;
                          else
                              return rhs.value < lhs.value;
                      });

            int sum = 0;
            int best = -2;

            for (int i = -1; i < r_size; ++i) {
                if (slv.bound_min(k) <= sum && sum <= slv.bound_max(k))
                    best = i;

                if (best != -2 && i + 1 < r_size &&
                    stop_iterating<Mode>(slv.R[i + 1].value))
                    break;

                auto var = begin + slv.R[i].id;
                sum += slv.factor(var->value);
            }

            int i = 0;
            for (; i <= best; ++i) {
                auto var = begin + slv.R[i].id;
                x.set(var->column);
            }

            for (; i != r_size; ++i) {
                auto var = begin + slv.R[i].id;
                x.unset(var->column);
            }
        }
    }

public:
    solver_initializer(Solver& slv,
                       solver_parameters::init_policy_type policy,
                       double init_policy_random,
                       double init_random)
      : dist(init_policy_random)
      , toss_up(init_random)
    {
        // x.clear();
        slv.reset();

        switch (policy) {
        case solver_parameters::init_policy_type::bastert:
            cycle = cycle_automaton::bastert;
            use_cycle = false;
            break;
        case solver_parameters::init_policy_type::pessimistic_solve:
            cycle = cycle_automaton::pessimistic_solve;
            use_cycle = false;
            break;
        case solver_parameters::init_policy_type::optimistic_solve:
            cycle = cycle_automaton::optimistic_solve;
            use_cycle = false;
            break;
        case solver_parameters::init_policy_type::cycle:
            cycle = cycle_automaton::pessimistic_solve;
            use_cycle = true;
            break;
        }
    }

    ~solver_initializer() noexcept
    {
        to_log(stdout,
               6u,
               "- solver-initializer status: bastert {} /"
               "pessimistic_solve {} / optimistic_solve {} / crossover {}\n",
               cycle_counter[0],
               cycle_counter[1],
               cycle_counter[2],
               cycle_counter[3]);
    }

    void next_state(bool is_a_solution) noexcept
    {
        if (is_a_solution) {
            single = single_automaton::init;
        } else {
            switch (single) {
            case single_automaton::init:
                single = single_automaton::improve_x_1;
                break;
            case single_automaton::improve_x_1:
                single = single_automaton::improve_x_2;
                break;
            case single_automaton::improve_x_2:
                single = single_automaton::improve_x_3;
                break;
            case single_automaton::improve_x_3:
                single = single_automaton::init;
            }
        }

        if (single == single_automaton::init) {
            if (use_cycle) {
                switch (cycle) {
                case cycle_automaton::bastert:
                    cycle = cycle_automaton::bastert_crossover;
                    break;
                case cycle_automaton::bastert_crossover:
                    cycle = cycle_automaton::pessimistic_solve;
                    break;
                case cycle_automaton::pessimistic_solve:
                    cycle = cycle_automaton::pessimistic_crossover;
                    break;
                case cycle_automaton::pessimistic_crossover:
                    cycle = cycle_automaton::optimistic_solve;
                    break;
                case cycle_automaton::optimistic_solve:
                    cycle = cycle_automaton::optimistic_crossover;
                    break;
                case cycle_automaton::optimistic_crossover:
                    cycle = cycle_automaton::bastert;
                    break;
                }
            } else {
                switch (cycle) {
                case cycle_automaton::bastert:
                    cycle = cycle_automaton::bastert_crossover;
                    break;
                case cycle_automaton::bastert_crossover:
                    cycle = cycle_automaton::bastert;
                    break;
                case cycle_automaton::pessimistic_solve:
                    cycle = cycle_automaton::pessimistic_crossover;
                    break;
                case cycle_automaton::pessimistic_crossover:
                    cycle = cycle_automaton::pessimistic_solve;
                    break;
                case cycle_automaton::optimistic_solve:
                    cycle = cycle_automaton::optimistic_crossover;
                    break;
                case cycle_automaton::optimistic_crossover:
                    cycle = cycle_automaton::optimistic_solve;
                    break;
                }
            }
        }
    }

    void init(Solver& slv, bit_array& x)
    {
        // x.clear();
        slv.reset();

        to_log(stdout,
               6u,
               "- initializer: init method {} in mode {} - init {}\n",
               cycle_automaton_string[static_cast<int>(cycle)],
               init_automaton_string[static_cast<int>(base)],
               single_automaton_string[static_cast<int>(single)]);

        for (int i = 0; i != slv.n; ++i)
            if (toss_up(slv.rng))
                x.set(i);
            else
                x.unset(i);

        switch (cycle) {
        case cycle_automaton::bastert:
            init_bastert(slv, x);
            break;
        case cycle_automaton::pessimistic_solve:
            init_pessimistic_solve(slv, x);
            break;
        case cycle_automaton::optimistic_solve:
            init_optimistic_solve(slv, x);
            break;
        case cycle_automaton::bastert_crossover:
        case cycle_automaton::pessimistic_crossover:
        case cycle_automaton::optimistic_crossover:
            bx_reach();
            break;
        }
    }

    Float reinit(Solver& slv,
                 bit_array& x,
                 bool is_a_solution,
                 const best_solution_recorder<Float, Mode>& best_recorder,
                 double kappa_min,
                 double kappa_max)
    {
        // x.clear();
        slv.reset();

        to_log(stdout,
               6u,
               "- initializer: re-init method {} in mode {} - init {} - "
               "is-a-solution {}\n",
               cycle_automaton_string[static_cast<int>(cycle)],
               init_automaton_string[static_cast<int>(base)],
               single_automaton_string[static_cast<int>(single)],
               is_a_solution);

        next_state(is_a_solution);

        switch (single) {
        case single_automaton::init:
            switch (base) {
            case init_automaton::random:
                for (int i = 0; i != slv.n; ++i)
                    if (toss_up(slv.rng))
                        x.set(i);
                    else
                        x.unset(i);

                base = init_automaton::current_best;
                break;
            case init_automaton::current_best:
                if (!best_recorder.copy_any_solution(x, slv.rng)) {
                    for (int i = 0; i != slv.n; ++i)
                        if (toss_up(slv.rng))
                            x.set(i);
                        else
                            x.unset(i);
                }

                base = init_automaton::best_recorder;
                break;
            case init_automaton::best_recorder:
                if (!best_recorder.copy_best_solution(x)) {
                    for (int i = 0; i != slv.n; ++i)
                        if (toss_up(slv.rng))
                            x.set(i);
                        else
                            x.unset(i);
                }

                base = init_automaton::random;
                break;
            }

            switch (cycle) {
            case cycle_automaton::bastert:
                init_bastert(slv, x);
                break;
            case cycle_automaton::pessimistic_solve:
                init_pessimistic_solve(slv, x);
                break;
            case cycle_automaton::optimistic_solve:
                init_optimistic_solve(slv, x);
                break;
            case cycle_automaton::bastert_crossover:
            case cycle_automaton::pessimistic_crossover:
            case cycle_automaton::optimistic_crossover:
                init_crossover(slv, x, best_recorder);
                break;
            }
            break;
        case single_automaton::improve_x_1:
            kappa_min += (kappa_max - kappa_min) / 20;
            break;
        case single_automaton::improve_x_2:
            kappa_min += (kappa_max - kappa_min) / 10;
            break;
        case single_automaton::improve_x_3:
            kappa_min += (kappa_max - kappa_min) / 5;
            break;
        }

        return static_cast<Float>(kappa_min);
    }
};

/**
 * Compute a problem lower or upper bounds based on Lagrangian multipliers
 * (valid if there are equality constraints only?)
 */
template<typename floatingpointT, typename modeT>
struct bounds_printer
{
    floatingpointT bestlb;
    floatingpointT bestub;
    floatingpointT max_cost;

    bounds_printer(floatingpointT max_cost_init)
      : bestlb(std::numeric_limits<floatingpointT>::lowest())
      , bestub(std::numeric_limits<floatingpointT>::max())
      , max_cost(max_cost_init)
    {}

    template<typename SolverT>
    floatingpointT init_bound(const SolverT& slv)
    {
        floatingpointT b{ 0 };

        for (auto c = 0; c != slv.m; ++c)
            b += slv.pi[c] * static_cast<floatingpointT>(slv.bound_init(c));

        return b;
    }

    template<typename SolverT>
    floatingpointT add_bound(const SolverT& slv,
                             int j,
                             floatingpointT sum_a_pi,
                             minimize_tag)
    {
        if (slv.c[j] - sum_a_pi < 0)
            return slv.c[j] - sum_a_pi;

        return { 0 };
    }

    template<typename SolverT>
    floatingpointT add_bound(const SolverT& slv,
                             int j,
                             floatingpointT sum_a_pi,
                             maximize_tag)
    {
        if (slv.c[j] - sum_a_pi > 0)
            return slv.c[j] - sum_a_pi;

        return { 0 };
    }

    void print_bound(const context_ptr& ctx,
                     floatingpointT lower_bound,
                     floatingpointT upper_bound,
                     minimize_tag)
    {
        bool better_gap = (lower_bound > bestlb || upper_bound < bestub);

        if (upper_bound < bestub)
            bestub = upper_bound;

        if (lower_bound > bestlb)
            bestlb = lower_bound;

        if (better_gap) {
            if (bestub == static_cast<floatingpointT>(0.0))
                info(ctx, "  - Lower bound: {}   (gap: 0%)\n", bestlb);
            else
                info(ctx,
                     "  - Lower bound: {}   (gap: {}%)\n",
                     bestlb,
                     static_cast<floatingpointT>(100.) * (bestub - bestlb) /
                       bestub);
        }
    }

    void print_bound(const context_ptr& ctx,
                     floatingpointT lower_bound,
                     floatingpointT upper_bound,
                     maximize_tag)
    {
        bool better_gap = (lower_bound > bestlb || upper_bound < bestub);

        if (upper_bound < bestub)
            bestub = upper_bound;

        if (lower_bound > bestlb)
            bestlb = lower_bound;

        if (better_gap) {
            if (bestlb == static_cast<floatingpointT>(0.0))
                info(ctx, "  - Upper bound: {}   (gap: 0%)\n", bestub);
            else
                info(ctx,
                     "  - Upper bound: {}   (gap: {}%)\n",
                     bestub,
                     static_cast<floatingpointT>(100.) * (bestlb - bestub) /
                       bestlb);
        }
    }

    floatingpointT init_ub(minimize_tag)
    {
        return std::numeric_limits<floatingpointT>::max();
    }

    floatingpointT init_ub(maximize_tag)
    {
        return std::numeric_limits<floatingpointT>::lowest();
    }

    template<typename SolverT>
    void operator()(const SolverT& slv,
                    const context_ptr& ctx,
                    const result& best)
    {
        floatingpointT lb = init_bound(slv);
        floatingpointT ub = init_ub(modeT());

        if (not best.solutions.empty())
            ub = static_cast<floatingpointT>(best.solutions.back().value);

        for (auto j = 0; j != slv.n; ++j)
            lb += add_bound(slv, j, slv.compute_sum_A_pi(j), modeT());

        lb *= max_cost; // restore original cost

        print_bound(ctx, lb, ub, modeT());
    }
};

struct compute_order
{
    std::vector<int> R;
    std::vector<std::pair<int, int>> m_order;
    solver_parameters::constraint_order order;
    bool use_cycle;

    compute_order(const solver_parameters::constraint_order order_,
                  const int variable_number)
      : R(static_cast<std::vector<int>::size_type>(variable_number))
      , m_order(static_cast<std::vector<int>::size_type>(variable_number))
      , order(order_ == solver_parameters::constraint_order::cycle
                ? solver_parameters::constraint_order::none
                : order_)
      , use_cycle(order_ == solver_parameters::constraint_order::cycle)
    {}

    static solver_parameters::constraint_order next_state(
      solver_parameters::constraint_order current) noexcept
    {
        using int_type = typename std::underlying_type<
          solver_parameters::constraint_order>::type;

        auto int_current = static_cast<int_type>(current);
        auto int_max =
          static_cast<int_type>(solver_parameters::constraint_order::cycle);

        ++int_current;

        return static_cast<solver_parameters::constraint_order>(
          int_current >= int_max ? 0 : int_current);
    }

    template<typename Solver, typename Xtype>
    void init(const Solver& s, const Xtype& x)
    {
        switch (order) {
        case solver_parameters::constraint_order::infeasibility_decr:
        case solver_parameters::constraint_order::infeasibility_incr:
            infeasibility_local_compute_violated_constraints(s, x);
            break;
        case solver_parameters::constraint_order::pi_sign_change:
            std::iota(R.begin(), R.end(), 0);