//
// Created by henri on 20.11.23.
//

#ifndef IDOL_KNAPSACKCOVER_H
#define IDOL_KNAPSACKCOVER_H

#include "idol/optimizers/mixed-integer-optimization/branch-and-bound/callbacks/BranchAndBoundCallbackFactory.h"
#include "idol/optimizers/mixed-integer-optimization/branch-and-bound/callbacks/BranchAndBoundCallback.h"
#include "idol/optimizers/mixed-integer-optimization/branch-and-bound/nodes/DefaultNodeInfo.h"

namespace idol::Cuts {
    template<class NodeInfoT> class KnapsackCover;
}

template<class NodeInfoT = idol::DefaultNodeInfo>
class idol::Cuts::KnapsackCover: public BranchAndBoundCallbackFactory<NodeInfoT> {
    std::optional<bool> m_lifting;
    std::optional<bool> m_apply_to_tree_nodes;
    std::optional<unsigned int> m_max_pass_root_node;
    std::optional<double> m_max_cuts_factor;
public:
    class Strategy;
    BranchAndBoundCallback<NodeInfoT> *operator()() override;
    BranchAndBoundCallbackFactory<NodeInfoT> *clone() const override;

    KnapsackCover& with_lifting(bool t_value);

    KnapsackCover& with_tree_node_cuts(bool t_value);

    KnapsackCover& with_max_pass_root_node(unsigned int t_value);

    KnapsackCover& with_max_cuts_factor(double t_value);
};

template<class NodeInfoT>
idol::Cuts::KnapsackCover<NodeInfoT> &idol::Cuts::KnapsackCover<NodeInfoT>::with_max_pass_root_node(unsigned int t_value) {

    if (m_max_pass_root_node.has_value()) {
        throw Exception("Maximum pass at root node has already been configured.");
    }

    m_max_pass_root_node = t_value;

    return *this;
}

template<class NodeInfoT>
idol::Cuts::KnapsackCover<NodeInfoT> &idol::Cuts::KnapsackCover<NodeInfoT>::with_max_cuts_factor(double t_value) {

    if (m_max_cuts_factor.has_value()) {
        throw Exception("Maximum cuts factor has already been configured.");
    }

    m_max_cuts_factor = t_value;

    return *this;
}

template<class NodeInfoT>
idol::Cuts::KnapsackCover<NodeInfoT> &idol::Cuts::KnapsackCover<NodeInfoT>::with_tree_node_cuts(bool t_value) {

    if (m_apply_to_tree_nodes.has_value()) {
        throw Exception("Tree node cuts have already been configured.");
    }

    m_apply_to_tree_nodes = t_value;

    return *this;
}

template<class NodeInfoT>
idol::Cuts::KnapsackCover<NodeInfoT> &idol::Cuts::KnapsackCover<NodeInfoT>::with_lifting(bool t_value) {

    if (m_lifting.has_value()) {
        throw Exception("Lifting has already been configured.");
    }

    m_lifting = t_value;

    return *this;
}

template<class NodeInfoT>
class idol::Cuts::KnapsackCover<NodeInfoT>::Strategy : public BranchAndBoundCallback<NodeInfoT> {

    struct KnapsackStructure {

        struct Item {
            Var variable;
            int weight;
            double current_value = -1.;
            std::optional<int> cut_coefficient;
            Item(const Var& t_var, int t_weight);

            Item(const Item&) = default;
            Item(Item&&) = default;

            Item& operator=(const Item&) = default;
            Item& operator=(Item&&) = default;
        };

        std::vector<Item> items;
        int capacity;
        KnapsackStructure(std::vector<Item> t_items, int t_capacity);
    };

    using ItemIterator = typename std::vector<typename KnapsackStructure::Item>::iterator;

    class SetOfItems {
        ItemIterator m_begin;
        ItemIterator m_end;
    public:
        SetOfItems(ItemIterator t_begin, ItemIterator t_end) : m_begin(t_begin), m_end(t_end) {}

        ItemIterator begin() const { return m_begin; }
        ItemIterator end() const { return m_end; }
        unsigned int size() const { return std::distance(m_begin, m_end); }

        static SetOfItems merge_consecutive(const SetOfItems& t_a, const SetOfItems& t_b);
    };

    std::list<KnapsackStructure> m_knapsack_structures;
    const bool m_lifting;
    const bool m_apply_to_tree_nodes;
    const unsigned int m_max_pass_root_node;
    const double m_max_cuts_factor;
    unsigned int m_max_cuts = 0;
    unsigned int m_n_cuts = 0;
    unsigned int m_n_separations = 0;
protected:
    // Detect structure
    void detect_knapsack_structure(const Ctr& t_ctr);
    void add_knapsack_structure(const Row& t_row, CtrType t_type);
    bool has_only_binary_variables(const Row& t_row);
    bool has_only_single_signed_coefficients(const Row& t_row);

    // Separation
    void separate_cut(const KnapsackStructure& t_knapsack_structure);
    std::tuple<SetOfItems, SetOfItems, SetOfItems> fix_variables_equal_to_one_or_zero(const SetOfItems& t_set_of_items);
    void set_current_values(const SetOfItems& t_set_of_items);
    int sum_of_weights(const SetOfItems& t_set_of_items);
    std::pair<SetOfItems, SetOfItems> compute_initial_cover(const SetOfItems& t_N_1, const SetOfItems& t_N_free, const SetOfItems& t_N_0, int t_capacity);
    std::pair<SetOfItems, SetOfItems> solve_knapsack_approximately(const SetOfItems& t_set_of_items, int t_capacity);
    std::pair<SetOfItems, SetOfItems> swap_items(const SetOfItems& t_a, const SetOfItems& t_b);
    std::pair<SetOfItems, SetOfItems> compute_minimal_cover(const SetOfItems& t_initial_cover, int t_capacity);
    std::pair<SetOfItems, SetOfItems> partition_minimal_cover(const SetOfItems& t_minimal_cover);
    std::pair<SetOfItems, SetOfItems> partition_remaining_items(const SetOfItems& t_remaining_items);
    void sort_by_non_decreasing_weights(const SetOfItems& t_set_of_items);
    void sort_by_non_increasing_weights(const SetOfItems& t_set_of_items);
    int sequential_up_and_down_lifting(const SetOfItems& t_C1, const SetOfItems& t_C2, const SetOfItems &t_F, const SetOfItems& t_R, int t_capacity);
    double compute_cut_activity(const SetOfItems& t_set_of_items, int t_right_hand_side);
    TempCtr create_cut(const SetOfItems& t_set_of_items, int t_right_hand_side);

    void debug(const std::string& t_name, const SetOfItems& t_set_of_items, bool t_with_values = false, bool t_with_cut = false);

    void initialize() override;

    void operator()(CallbackEvent t_event) override;

    void log_after_termination() override;

public:
    Strategy(bool t_use_lifting,
             bool t_apply_to_tree_nodes,
             unsigned int t_max_pass_root_node,
             double t_max_cuts_factor);
};

template<class NodeInfoT>
void idol::Cuts::KnapsackCover<NodeInfoT>::Strategy::log_after_termination() {

    BranchAndBoundCallback<NodeInfoT>::log_after_termination();

    std::cout << "\tKnapsack Cover Cuts: " << m_n_cuts << std::endl;
}

template<class NodeInfoT>
idol::Cuts::KnapsackCover<NodeInfoT>::Strategy::Strategy(bool t_use_lifting,
                                                         bool t_apply_to_tree_nodes,
                                                         unsigned int t_max_pass_root_node,
                                                         double t_max_cuts_factor)
    : m_lifting(t_use_lifting),
      m_apply_to_tree_nodes(t_apply_to_tree_nodes),
      m_max_pass_root_node(t_max_pass_root_node),
      m_max_cuts_factor(t_max_cuts_factor) {

}

template<class NodeInfoT>
void idol::Cuts::KnapsackCover<NodeInfoT>::Strategy::debug(const std::string &t_name,
                                                           const SetOfItems &t_set_of_items,
                                                           bool t_with_values,
                                                           bool t_with_cut) {

    std::cout << t_name << ": ";
    for (const auto& item : t_set_of_items) {
        if (t_with_cut) {
            if (item.cut_coefficient.has_value()) {
                std::cout << item.cut_coefficient.value();
            } else {
                std::cout << "?";
            }
            std::cout << " ";
        }
        std::cout << item.variable;
        if (t_with_values) {
            std::cout << "(" << item.current_value << ")";
        }
        std::cout << ", ";
    }
    std::cout << std::endl;

}

template<class NodeInfoT>
typename idol::Cuts::KnapsackCover<NodeInfoT>::Strategy::SetOfItems
idol::Cuts::KnapsackCover<NodeInfoT>::Strategy::SetOfItems::merge_consecutive(const SetOfItems &t_a,
                                                                             const SetOfItems &t_b) {

    if (t_a.end() != t_b.begin()) {
        throw Exception("Trying to merge non-consecutive sets of items while calling merge_consecutive.");
    }

    return {
        t_a.begin(),
        t_b.end()
    };

}

template<class NodeInfoT>
idol::Cuts::KnapsackCover<NodeInfoT>::Strategy::KnapsackStructure::KnapsackStructure(std::vector<Item> t_items, int t_capacity)
        : items(std::move(t_items)), capacity(t_capacity) {

}

template<class NodeInfoT>
idol::Cuts::KnapsackCover<NodeInfoT>::Strategy::KnapsackStructure::Item::Item(const idol::Var &t_var, int t_weight)
        : variable(t_var), weight(t_weight) {

}

template<class NodeInfoT>
void idol::Cuts::KnapsackCover<NodeInfoT>::Strategy::initialize() {
    BranchAndBoundCallback<NodeInfoT>::initialize();

    m_knapsack_structures.clear();
    m_n_cuts = 0;
    m_n_separations = 0;

    const auto& model = this->original_model();

    for (const auto& ctr : model.ctrs()) {
        detect_knapsack_structure(ctr);
    }

    m_max_cuts = m_max_cuts_factor * model.ctrs().size();

    std::cout << "Lifted Minimal Cover Cuts 1, " << m_knapsack_structures.size() << " candidates found" << std::endl;
}

template<class NodeInfoT>
void idol::Cuts::KnapsackCover<NodeInfoT>::Strategy::detect_knapsack_structure(const idol::Ctr &t_ctr) {

    const auto& model = this->original_model();
    const auto& row = model.get_ctr_row(t_ctr);

    if (!row.quadratic().empty()) {
        return;
    }

    if (!has_only_binary_variables(row)) {
        return;
    }

    if (!has_only_single_signed_coefficients(row)) {
        return;
    }

    const auto type = model.get_ctr_type(t_ctr);

    if (type == Equal) {
        return;
    }

    add_knapsack_structure(row, type);

}

template<class NodeInfoT>
void idol::Cuts::KnapsackCover<NodeInfoT>::Strategy::add_knapsack_structure(const idol::Row &t_row, idol::CtrType t_type) {

    const unsigned int n_items = t_row.linear().size();
    const double factor = (t_type == LessOrEqual ? 1. : -1.);

    auto copy = t_row;
    copy.scale_to_integers(Tolerance::Digits);

    std::vector<typename KnapsackStructure::Item> items;
    items.reserve(n_items);

    for (const auto& [var, coefficient] : copy.linear()) {
        items.emplace_back(var, factor * coefficient.as_numerical());
    }

    m_knapsack_structures.emplace_back(std::move(items), factor * copy.rhs().as_numerical());

}

template<class NodeInfoT>
bool idol::Cuts::KnapsackCover<NodeInfoT>::Strategy::has_only_binary_variables(const idol::Row &t_row) {

    const auto& model = this->original_model();

    for (const auto& [var, coefficient] : t_row.linear()) {

        if (model.get_var_type(var) != Binary) {
            return false;
        }

    }

    return true;
}

template<class NodeInfoT>
bool idol::Cuts::KnapsackCover<NodeInfoT>::Strategy::has_only_single_signed_coefficients(const idol::Row &t_row) {

    std::optional<bool> are_all_non_negative;

    for (const auto& [var, coefficient] : t_row.linear()) {

        bool is_non_negative = coefficient.as_numerical() >= 0;

        if (!are_all_non_negative.has_value()) [[unlikely]] {
            are_all_non_negative = is_non_negative;
            continue;
        }

        if (is_non_negative != are_all_non_negative) {
            return false;
        }

    }

    return true;
}

template<class NodeInfoT>
void idol::Cuts::KnapsackCover<NodeInfoT>::Strategy::operator()(idol::CallbackEvent t_event) {

    if (t_event != InvalidSolution) {
        return;
    }

    if (m_n_cuts >= m_max_cuts) {
        return;
    }

    if (this->node().level() == 0) {

        if (m_n_separations > m_max_pass_root_node) {
            return;
        }

    }

    if (!m_apply_to_tree_nodes && this->node().level() != 0) {
        return;
    }

    for (const auto& knapsack_structure : m_knapsack_structures) {
        separate_cut(knapsack_structure);
        ++m_n_separations;
    }

}

template<class NodeInfoT>
void idol::Cuts::KnapsackCover<NodeInfoT>::Strategy::separate_cut(const idol::Cuts::KnapsackCover<NodeInfoT>::Strategy::KnapsackStructure &t_knapsack_structure) {

    if (m_n_cuts >= m_max_cuts) {
        return;
    }

    auto knapsack_structure = t_knapsack_structure;
    SetOfItems N(knapsack_structure.items.begin(), knapsack_structure.items.end());

    set_current_values(N);

    const auto [N_1, N_free, N_0] = fix_variables_equal_to_one_or_zero(N);
    const int sum_weights_of_variables_not_fixed_to_zero = sum_of_weights(N_1) + sum_of_weights(N_free);

    if (sum_weights_of_variables_not_fixed_to_zero < knapsack_structure.capacity + 1) {
        return;
    }

    // Initial cover
    auto [initial_cover, not_initial_cover] = compute_initial_cover(
                                                                N_1,
                                                                N_free,
                                                                N_0,
                                                                sum_weights_of_variables_not_fixed_to_zero - (knapsack_structure.capacity + 1)
                                                                );

    // Minimal cover
    auto [C, others] = compute_minimal_cover(initial_cover, t_knapsack_structure.capacity);

    int rhs;

    if (m_lifting) {

        // Partition
        auto not_in_C = SetOfItems::merge_consecutive(others, not_initial_cover);
        auto [C1, C2] = partition_minimal_cover(C);

        // Lifting sequence and computing the lifting coefficients
        auto [F, R] = partition_remaining_items(not_in_C);
        sort_by_non_decreasing_weights(C1);
        sort_by_non_increasing_weights(C2);
        sort_by_non_increasing_weights(F);
        sort_by_non_increasing_weights(R);

        rhs = sequential_up_and_down_lifting(C1, C2, F, R, knapsack_structure.capacity);

    } else {

        rhs = C.size() - 1;
        for (auto& item : C) {
            item.cut_coefficient = 1.;
        }

    }

    const double activity = compute_cut_activity(N, rhs);

    if (activity < Tolerance::Feasibility) {
        return;
    }

    auto cut = create_cut(N, rhs);

    this->add_user_cut(cut);
    ++m_n_cuts;

}

template<class NodeInfoT>
idol::TempCtr idol::Cuts::KnapsackCover<NodeInfoT>::Strategy::create_cut(const SetOfItems &t_set_of_items, int t_right_hand_side) {

    Expr result;

    for (const auto& item : t_set_of_items) {

        if (!item.cut_coefficient.has_value()) {
            continue;
        }

        result += item.cut_coefficient.value() * item.variable;

    }

    return result <= t_right_hand_side;
}

template<class NodeInfoT>
double idol::Cuts::KnapsackCover<NodeInfoT>::Strategy::compute_cut_activity(const SetOfItems &t_set_of_items, int t_right_hand_side) {

    double result = 0.;

    for (const auto& item : t_set_of_items) {

        if (!item.cut_coefficient.has_value()) {
            continue;
        }

        result += item.cut_coefficient.value() * item.current_value;

    }

    return result - t_right_hand_side;
}

template<class NodeInfoT>
int idol::Cuts::KnapsackCover<NodeInfoT>::Strategy::sequential_up_and_down_lifting(const SetOfItems &t_C1,
                                                                                    const SetOfItems &t_C2,
                                                                                    const SetOfItems &t_F,
                                                                                    const SetOfItems &t_R,
                                                                                    int t_capacity) {

    // Set coefficients of variables in C1 to one
    for (auto& item : t_C1) {
        item.cut_coefficient = 1;
    }

    const auto initialize_min_weights = [](const SetOfItems& t_C1) {

        const unsigned int n_C1 = t_C1.size();

        std::vector<int> result;
        result.reserve(n_C1 + 1);
        result.emplace_back(0);

        auto it = t_C1.begin();
        for (unsigned int w = 0 ; w < n_C1 ; ++w) {
            result.emplace_back(result[w] + it->weight);
            ++it;
        }

        return result;
    };

    auto min_weights = initialize_min_weights(t_C1);

    int alpha_0 = t_C1.size() - 1;
    int lift_rhs = alpha_0;

    int fixed_ones_weight = sum_of_weights(t_C2);

    assert(fixed_ones_weight >= 0);

    // Lifting items in F
    for (auto& item : t_F) {

        int z;
        if (t_capacity - fixed_ones_weight - item.weight < 0) { // knapsack is infeasible
            z = 0;
        } else if (min_weights[alpha_0] <= t_capacity - fixed_ones_weight - item.weight) { // easy case
            z = lift_rhs;
        } else { // binary search for z

            int lb = 0;
            int ub = lift_rhs + 1;
            while (lb < ub - 1) {
                const int middle = (lb + ub) / 2;
                if ( min_weights[middle] <= t_capacity - fixed_ones_weight - item.weight ) {
                    lb = middle;
                } else {
                    ub = middle;
                }
            }

            assert(lb == ub - 1);
            assert(0 <= lb && lb < min_weights.size());
            assert(min_weights[lb] <= t_capacity - fixed_ones_weight - item.weight);
            assert(lb == min_weights.size() - 1 || min_weights[lb + 1] > t_capacity - fixed_ones_weight - item.weight);

            z = lb;

            assert(z <= lift_rhs);

        }

        const int alpha_j = alpha_0 - z;

        item.cut_coefficient = alpha_j;

        if (alpha_j == 0) {
            continue;
        }

        min_weights.resize(min_weights.size() + alpha_j, std::numeric_limits<int>::infinity());

        for (int w = min_weights.size() - 1 ; w >= 0 ; w--) {

            if (w < alpha_j ) {
                min_weights[w] = std::min(min_weights[w], item.weight);
            } else {
                min_weights[w] = std::min(min_weights[w], min_weights[w - alpha_j] + item.weight);
            }

        }

    }

    // Lifting in C2
    for (auto& item : t_C2) {

        int lb = 0;
        int ub = min_weights.size();
        while (lb < ub - 1) {
            int middle = (lb + ub) / 2;
            if ( min_weights[middle] <= t_capacity - fixed_ones_weight + item.weight ) {
                lb = middle;
            } else {
                ub = middle;
            }
        }
        int z = lb;

        const int alpha_j = z - alpha_0;

        item.cut_coefficient = alpha_j;

        alpha_0 += alpha_j;

        fixed_ones_weight -= item.weight;

        if (alpha_j == 0) {
            continue;
        }

        if (alpha_j > 0) {
            min_weights.resize(min_weights.size() + alpha_j, std::numeric_limits<int>::infinity());
        }

        for (int w = min_weights.size() - 1 ; w >= 0 ; w--) {

            if (w < alpha_j ) {
                min_weights[w] = std::min(min_weights[w], item.weight);
            } else {
                min_weights[w] = std::min(min_weights[w], min_weights[w - alpha_j] + item.weight);
            }

        }


    }

    assert(fixed_ones_weight == 0);

    // Lifting in R
    for (auto& item : t_R) {

        int z;
        if (min_weights[alpha_0] <= t_capacity - item.weight) {
            z = alpha_0;
        } else { // binary search for z

            int lb = 0;
            int ub = alpha_0 + 1;
            while (lb < ub - 1) {
                int middle = (lb + ub) / 2;
                if ( min_weights[middle] <= t_capacity - item.weight ) {
                    lb = middle;
                } else {
                    ub = middle;
                }
            }
            z = lb;

        }

        const int alpha_j = alpha_0 - z;

        item.cut_coefficient = alpha_j;

        if (alpha_j == 0) {
            continue;
        }

        for (int w = min_weights.size() - 1 ; w >= 0 ; w--) {

            if (w < alpha_j ) {
                min_weights[w] = std::min(min_weights[w], item.weight);
            } else {
                min_weights[w] = std::min(min_weights[w], min_weights[w - alpha_j] + item.weight);
            }

        }


    }

    return alpha_0;

}

template<class NodeInfoT>
void idol::Cuts::KnapsackCover<NodeInfoT>::Strategy::sort_by_non_increasing_weights(const SetOfItems &t_set_of_items) {

    std::sort(t_set_of_items.begin(), t_set_of_items.end(), [](const auto& t_a, const auto& t_b) {
        return t_a.weight >= t_b.weight;
    });

}

template<class NodeInfoT>
void idol::Cuts::KnapsackCover<NodeInfoT>::Strategy::sort_by_non_decreasing_weights(const SetOfItems &t_set_of_items) {

    std::sort(t_set_of_items.begin(), t_set_of_items.end(), [](const auto& t_a, const auto& t_b) {
        return t_a.weight <= t_b.weight;
    });

}

template<class NodeInfoT>
void idol::Cuts::KnapsackCover<NodeInfoT>::Strategy::set_current_values(const SetOfItems &t_set_of_items) {

    const auto& primal_solution = this->node().info().primal_solution();

    for (auto& item : t_set_of_items) {
        item.current_value = primal_solution.get(item.variable);
    }

}

template<class NodeInfoT>
std::tuple<
        typename idol::Cuts::KnapsackCover<NodeInfoT>::Strategy::SetOfItems,
        typename idol::Cuts::KnapsackCover<NodeInfoT>::Strategy::SetOfItems,
        typename idol::Cuts::KnapsackCover<NodeInfoT>::Strategy::SetOfItems
        >
idol::Cuts::KnapsackCover<NodeInfoT>::Strategy::fix_variables_equal_to_one_or_zero(const SetOfItems &t_set_of_items) {

    auto end_of_fixed_to_one = t_set_of_items.begin();
    auto begin_of_fixed_to_zero = t_set_of_items.end();

    for (auto it = end_of_fixed_to_one ; it != begin_of_fixed_to_zero && end_of_fixed_to_one != begin_of_fixed_to_zero ; ) {

        if (equals(it->current_value, 0., Tolerance::Integer)) {
            --begin_of_fixed_to_zero;
            std::iter_swap(it, begin_of_fixed_to_zero);
            continue;
        }

        if (equals(it->current_value, 1., Tolerance::Integer)) {
            std::iter_swap(it, end_of_fixed_to_one);
            ++end_of_fixed_to_one;
            ++it;
            continue;
        }

        ++it;

    }

    return {
            { t_set_of_items.begin(), end_of_fixed_to_one }, // fixed to 1
            { end_of_fixed_to_one, begin_of_fixed_to_zero }, // free
            { begin_of_fixed_to_zero, t_set_of_items.end() }, // fixed to 0
    };
}

template<class NodeInfoT>
int idol::Cuts::KnapsackCover<NodeInfoT>::Strategy::sum_of_weights(const SetOfItems &t_set_of_items) {

    int result = 0.;

    for (auto it = t_set_of_items.begin(), end = t_set_of_items.end() ; it != end ; ++it) {
        result += it->weight;
    }

    return result;
}

template<class NodeInfoT>
std::pair<
        typename idol::Cuts::KnapsackCover<NodeInfoT>::Strategy::SetOfItems,
        typename idol::Cuts::KnapsackCover<NodeInfoT>::Strategy::SetOfItems
>
idol::Cuts::KnapsackCover<NodeInfoT>::Strategy::compute_initial_cover(const SetOfItems &t_N_1,
                                                                     const SetOfItems &t_N_free,
                                                                     const SetOfItems &t_N_0,
                                                                     int t_capacity) {

    auto [ones, zeroes] = solve_knapsack_approximately(t_N_free, t_capacity);

    assert(ones.end() == zeroes.begin());

    auto [flipped_zeroes, flipped_ones] = swap_items(ones, zeroes);

    auto initial_cover = SetOfItems::merge_consecutive(t_N_1, flipped_zeroes);
    auto not_initial_cover = SetOfItems::merge_consecutive(flipped_ones, t_N_0);

    return { initial_cover, not_initial_cover };

}

template<class NodeInfoT>
std::pair<
        typename idol::Cuts::KnapsackCover<NodeInfoT>::Strategy::SetOfItems,
        typename idol::Cuts::KnapsackCover<NodeInfoT>::Strategy::SetOfItems
>
idol::Cuts::KnapsackCover<NodeInfoT>::Strategy::partition_minimal_cover(const SetOfItems &t_minimal_cover) {

    auto end_of_C1 = t_minimal_cover.end();

    for (auto it = t_minimal_cover.begin() ; it != end_of_C1 ; ) {

        if (!equals(it->current_value, 1., Tolerance::Integer)) {
            ++it;
            continue;
        }

        --end_of_C1;
        std::iter_swap(it, end_of_C1);

    }

    return {
            { t_minimal_cover.begin(), end_of_C1 },
            { end_of_C1, t_minimal_cover.end() }
    };
}

template<class NodeInfoT>
std::pair<
        typename idol::Cuts::KnapsackCover<NodeInfoT>::Strategy::SetOfItems,
        typename idol::Cuts::KnapsackCover<NodeInfoT>::Strategy::SetOfItems
>
idol::Cuts::KnapsackCover<NodeInfoT>::Strategy::partition_remaining_items(const SetOfItems &t_remaining_items) {

    auto end_of_F = t_remaining_items.end();

    for (auto it = t_remaining_items.begin() ; it != end_of_F ; ) {

        if (!equals(it->current_value, 0, Tolerance::Integer)) {
            ++it;
            continue;
        }

        --end_of_F;
        std::iter_swap(it, end_of_F);

    }

    return {
            { t_remaining_items.begin(), end_of_F },
            { end_of_F, t_remaining_items.end() },
    };
}

template<class NodeInfoT>
std::pair<
        typename idol::Cuts::KnapsackCover<NodeInfoT>::Strategy::SetOfItems,
        typename idol::Cuts::KnapsackCover<NodeInfoT>::Strategy::SetOfItems
    >
idol::Cuts::KnapsackCover<NodeInfoT>::Strategy::solve_knapsack_approximately(const SetOfItems &t_set_of_items,
                                                                            int t_capacity) {

    std::sort(t_set_of_items.begin(),
              t_set_of_items.end(),
              [&](const auto& t_a, const auto& t_b) {

        const double score_a = (1. - t_a.current_value) / t_a.weight;
        const double score_b = (1. - t_b.current_value) / t_b.weight;

        return score_a >= score_b;

    });

    int weight_of_packed_items = 0.;
    auto it = t_set_of_items.begin();
    auto end = t_set_of_items.end();

    auto begin_zero = it;

    for ( ; it != end ; ++it) {

        if (weight_of_packed_items + it->weight > t_capacity) {
            // all remaining-s go to zero
            break;
        }

        // set to 1
        std::iter_swap(it, begin_zero);
        ++begin_zero;

        weight_of_packed_items += it->weight;

    }

    return {
            { t_set_of_items.begin(), begin_zero }, // ones
            { begin_zero, t_set_of_items.end() }, // zeroes
    };

}

template<class NodeInfoT>
std::pair<
        typename idol::Cuts::KnapsackCover<NodeInfoT>::Strategy::SetOfItems,
        typename idol::Cuts::KnapsackCover<NodeInfoT>::Strategy::SetOfItems
>
idol::Cuts::KnapsackCover<NodeInfoT>::Strategy::swap_items(const SetOfItems &t_a,
                                                          const SetOfItems &t_b) {

    const auto min_size = std::min(t_a.size(), t_b.size());
    const auto max_size = std::max(t_a.size(), t_b.size());

    std::swap_ranges(t_a.begin(), t_a.begin() + min_size, t_a.begin() + max_size);

    ItemIterator split = t_a.begin() + t_b.size();

    return {
            { t_a.begin(), split },
            { split, t_b.end() }
    };
}

template<class NodeInfoT>
std::pair<
    typename idol::Cuts::KnapsackCover<NodeInfoT>::Strategy::SetOfItems,
    typename idol::Cuts::KnapsackCover<NodeInfoT>::Strategy::SetOfItems
>
idol::Cuts::KnapsackCover<NodeInfoT>::Strategy::compute_minimal_cover(const SetOfItems &t_initial_cover,
                                                                     int t_capacity) {

    std::sort(t_initial_cover.begin(), t_initial_cover.end(), [](const auto& t_a, const auto& t_b) {

        const double score_a = (1. - t_a.current_value) / t_a.weight;
        const double score_b = (1. - t_b.current_value) / t_b.weight;

        if (equals(score_a, score_b, Tolerance::Sparsity)) {
            return t_a.weight <= t_b.weight;
        }

        return score_a >= score_b;
    });

    int sum_of_weights_in_C = sum_of_weights(t_initial_cover);
    auto end = t_initial_cover.end();

    for (auto it = t_initial_cover.begin() ; it != end ; ) {

        if (sum_of_weights_in_C - it->weight <= t_capacity) {
            ++it;
            continue;
        }

        sum_of_weights_in_C -= it->weight;
        --end;
        std::iter_swap(it, end);

    }

    return {
        { t_initial_cover.begin(), end },
        { end, t_initial_cover.end() },
    };

}

template<class NodeInfoT>
idol::BranchAndBoundCallbackFactory<NodeInfoT> *idol::Cuts::KnapsackCover<NodeInfoT>::clone() const {
    return new KnapsackCover<NodeInfoT>(*this);
}

template<class NodeInfoT>
idol::BranchAndBoundCallback<NodeInfoT> *idol::Cuts::KnapsackCover<NodeInfoT>::operator()() {
    return new Strategy(
            m_lifting.has_value() ? m_lifting.value() : true,
            m_apply_to_tree_nodes.has_value() ? m_apply_to_tree_nodes.value() : true,
            m_max_pass_root_node.has_value() ? m_max_pass_root_node.value() : 200,
            m_max_cuts_factor.has_value() ? m_max_cuts_factor.value() : 1.5
        );
}

#endif //IDOL_KNAPSACKCOVER_H
