Integer Linear Programming
(optimization/integer_lp.hpp)

View this file on GitHub
Last update: 2026-07-07 14:26:59+09:00
Include: #include "optimization/integer_lp.hpp"

Overview

integer_lp_maximize(a, b, c) solves an integer linear programming problem in standard inequality form:

\[\begin{array}{ll} \text{maximize} & c^T x \\ \text{subject to} & A x \le b \\ & x \in \mathbb{Z}_{\ge 0}^W \end{array}\]

The implementation uses branch-and-bound over LP relaxations solved by simplex_maximize. It is useful for small or naturally bounded instances. Integer linear programming is NP-hard, so the number of explored nodes can be exponential.

integer_lp_minimize(a, b, c) solves the corresponding minimization problem with the same constraints. integer_lp(a, b, c) is an alias of integer_lp_maximize(a, b, c).

Interface

The inputs are:

Argument	Type	Meaning
`a`	`std::vector<std::vector<T>>`	Constraint matrix `A`.
`b`	`std::vector<T>`	Right-hand side vector. Constraint `i` is `a[i] * x <= b[i]`.
`c`	`std::vector<T>`	Objective coefficients.
`eps`	`long double`	Optional tolerance used by LP relaxations. The default is `1e-10`.

a.size() must equal b.size(), and every row of a must have c.size() entries. T must be a signed integer type, such as int or long long.

IntegerLpStatus has these values:

Value	Meaning
`IntegerLpStatus::Optimal`	A finite integer optimum was found.
`IntegerLpStatus::Infeasible`	No integer vector satisfies all constraints.
`IntegerLpStatus::Unbounded`	The objective is unbounded in the requested direction.

IntegerLpResult<T> contains these members:

Member / Method	Type / Signature	Meaning
`status`	`IntegerLpStatus`	Solver status.
`objective_value`	`T`	Optimal objective value when `status` is `Optimal`.
`variables`	`std::vector<T>`	Optimal variable values when `status` is `Optimal`.
`is_optimal`	`bool is_optimal() const`	Returns whether the status is `Optimal`.
`is_infeasible`	`bool is_infeasible() const`	Returns whether the status is `Infeasible`.
`is_unbounded`	`bool is_unbounded() const`	Returns whether the status is `Unbounded`.

When the result is Unbounded, variables contains an integer feasible point found during the search, and objective_value is set to the numeric limit in the unbounded direction. For Infeasible, objective_value and variables are placeholders.

Functions

Function	Signature	Description	Complexity
`integer_lp_maximize`	`template <class T> IntegerLpResult<T> integer_lp_maximize(const std::vector<std::vector<T>>& a, const std::vector<T>& b, const std::vector<T>& c, long double eps = 1e-10L)`	Maximizes `c^T x` subject to `A x <= b` and nonnegative integer `x`.	Exponential in the worst case.
`integer_lp_minimize`	`template <class T> IntegerLpResult<T> integer_lp_minimize(const std::vector<std::vector<T>>& a, const std::vector<T>& b, const std::vector<T>& c, long double eps = 1e-10L)`	Minimizes `c^T x` under the same constraints.	Exponential in the worst case.
`integer_lp`	`template <class T> IntegerLpResult<T> integer_lp(const std::vector<std::vector<T>>& a, const std::vector<T>& b, const std::vector<T>& c, long double eps = 1e-10L)`	Alias of `integer_lp_maximize`.	Exponential in the worst case.

Example

#include "optimization/integer_lp.hpp"
#include <iostream>
#include <vector>

int main() {
    std::vector<std::vector<long long>> a;
    a.emplace_back(std::vector<long long>{2, 1});
    a.emplace_back(std::vector<long long>{1, 2});

    std::vector<long long> b = {4, 4};
    std::vector<long long> c = {3, 2};

    auto result = m1une::opt::integer_lp_maximize(a, b, c);
    if (result.is_optimal()) {
        std::cout << result.objective_value << "\n";  // 6
        std::cout << result.variables[0] << " " << result.variables[1] << "\n";
    }
}

Depends on

Simplex Algorithm (optimization/simplex.hpp)

Required by

Optimization All (optimization/all.hpp)

Verified with

Code

#ifndef M1UNE_OPTIMIZATION_INTEGER_LP_HPP
#define M1UNE_OPTIMIZATION_INTEGER_LP_HPP 1

#include <algorithm>
#include <cassert>
#include <cmath>
#include <limits>
#include <type_traits>
#include <vector>

#include "simplex.hpp"

namespace m1une {
namespace opt {

enum class IntegerLpStatus {
    Optimal,
    Infeasible,
    Unbounded,
};

template <class T>
struct IntegerLpResult {
    IntegerLpStatus status;
    T objective_value;
    std::vector<T> variables;

    bool is_optimal() const { return status == IntegerLpStatus::Optimal; }
    bool is_infeasible() const { return status == IntegerLpStatus::Infeasible; }
    bool is_unbounded() const { return status == IntegerLpStatus::Unbounded; }
};

namespace detail {

template <class T>
struct IntegerLpSolver {
    using Real = long double;

    struct Node {
        std::vector<std::vector<Real>> a;
        std::vector<Real> b;
    };

    int variable_count;
    bool maximize;
    Real eps;
    std::vector<T> objective;
    std::vector<Real> relaxation_objective;
    Node initial_node;

    bool has_incumbent = false;
    T best_value = T();
    std::vector<T> best_variables;

    IntegerLpSolver(const std::vector<std::vector<T>>& a, const std::vector<T>& b,
                    const std::vector<T>& c, bool is_maximize, Real epsilon)
        : variable_count(int(c.size())),
          maximize(is_maximize),
          eps(epsilon),
          objective(c),
          relaxation_objective(c.size(), Real()),
          initial_node() {
        initial_node.a.assign(a.size(), std::vector<Real>(variable_count, Real()));
        initial_node.b.assign(b.size(), Real());
        for (int i = 0; i < int(a.size()); i++) {
            for (int j = 0; j < variable_count; j++) initial_node.a[i][j] = Real(a[i][j]);
            initial_node.b[i] = Real(b[i]);
        }
        Real sign = maximize ? Real(1) : Real(-1);
        for (int j = 0; j < variable_count; j++) relaxation_objective[j] = sign * Real(c[j]);
    }

    Real abs_value(Real x) const {
        return x < Real() ? -x : x;
    }

    bool better_value(T lhs, T rhs) const {
        return maximize ? lhs > rhs : lhs < rhs;
    }

    bool can_prune_by_bound(Real relaxation_value) const {
        if (!has_incumbent) return false;
        Real signed_best = maximize ? Real(best_value) : -Real(best_value);
        return relaxation_value <= signed_best + eps;
    }

    T evaluate(const std::vector<T>& variables) const {
        T result = T();
        for (int i = 0; i < variable_count; i++) result += objective[i] * variables[i];
        return result;
    }

    bool round_solution(const std::vector<Real>& real_variables, std::vector<T>& variables) const {
        variables.assign(variable_count, T());
        for (int i = 0; i < variable_count; i++) {
            Real value = real_variables[i];
            if (value < -eps) return false;
            Real rounded = std::round(value);
            if (abs_value(value - rounded) > eps) return false;
            variables[i] = static_cast<T>(rounded);
        }
        return true;
    }

    int find_fractional_variable(const std::vector<Real>& real_variables) const {
        int result = -1;
        Real best_distance = eps;
        for (int i = 0; i < variable_count; i++) {
            Real value = real_variables[i];
            Real rounded = std::round(value);
            Real distance = abs_value(value - rounded);
            if (distance > best_distance) {
                best_distance = distance;
                result = i;
            }
        }
        return result;
    }

    Node with_upper_bound(const Node& node, int variable, T bound) const {
        Node result = node;
        result.a.emplace_back(variable_count, Real());
        result.a.back()[variable] = Real(1);
        result.b.push_back(Real(bound));
        return result;
    }

    Node with_lower_bound(const Node& node, int variable, T bound) const {
        Node result = node;
        result.a.emplace_back(variable_count, Real());
        result.a.back()[variable] = Real(-1);
        result.b.push_back(-Real(bound));
        return result;
    }

    void push_branches(std::vector<Node>& stack, const Node& node, int variable, Real value) const {
        Real floor_value = std::floor(value);
        Real ceil_value = std::ceil(value);
        T upper_bound = static_cast<T>(floor_value);
        T lower_bound = static_cast<T>(ceil_value);

        bool has_upper_branch = upper_bound >= T();
        bool prefer_lower_branch = relaxation_objective[variable] >= -eps;

        if (prefer_lower_branch) {
            if (has_upper_branch) stack.push_back(with_upper_bound(node, variable, upper_bound));
            stack.push_back(with_lower_bound(node, variable, lower_bound));
        } else {
            stack.push_back(with_lower_bound(node, variable, lower_bound));
            if (has_upper_branch) stack.push_back(with_upper_bound(node, variable, upper_bound));
        }
    }

    bool has_positive_direction(const Node& node) const {
        std::vector<std::vector<Real>> direction_a = node.a;
        std::vector<Real> direction_b(node.b.size(), Real());

        std::vector<Real> objective_row(variable_count, Real());
        for (int i = 0; i < variable_count; i++) objective_row[i] = -relaxation_objective[i];
        direction_a.push_back(objective_row);
        direction_b.push_back(Real(-1));

        std::vector<Real> zero_objective(variable_count, Real());
        auto result = simplex_maximize(direction_a, direction_b, zero_objective, eps);
        return result.is_optimal();
    }

    bool find_integer_feasible(const Node& start, std::vector<T>& feasible_variables) const {
        std::vector<Node> stack;
        stack.push_back(start);
        std::vector<Real> zero_objective(variable_count, Real());

        while (!stack.empty()) {
            Node node = stack.back();
            stack.pop_back();

            auto relaxation = simplex_maximize(node.a, node.b, zero_objective, eps);
            if (relaxation.is_infeasible()) continue;
            if (relaxation.is_unbounded()) continue;

            if (round_solution(relaxation.variables, feasible_variables)) return true;

            int variable = find_fractional_variable(relaxation.variables);
            if (variable == -1) continue;
            push_branches(stack, node, variable, relaxation.variables[variable]);
        }
        return false;
    }

    void update_incumbent(const std::vector<T>& variables) {
        T value = evaluate(variables);
        if (!has_incumbent || better_value(value, best_value)) {
            has_incumbent = true;
            best_value = value;
            best_variables = variables;
        }
    }

    IntegerLpResult<T> make_infeasible_result() const {
        IntegerLpResult<T> result;
        result.status = IntegerLpStatus::Infeasible;
        result.objective_value = T();
        result.variables.assign(variable_count, T());
        return result;
    }

    IntegerLpResult<T> make_unbounded_result(const std::vector<T>& variables) const {
        IntegerLpResult<T> result;
        result.status = IntegerLpStatus::Unbounded;
        result.objective_value =
            maximize ? std::numeric_limits<T>::max() : std::numeric_limits<T>::lowest();
        result.variables = variables;
        return result;
    }

    IntegerLpResult<T> make_optimal_result() const {
        IntegerLpResult<T> result;
        result.status = IntegerLpStatus::Optimal;
        result.objective_value = best_value;
        result.variables = best_variables;
        return result;
    }

    IntegerLpResult<T> solve() {
        std::vector<Node> stack;
        stack.push_back(initial_node);

        while (!stack.empty()) {
            Node node = stack.back();
            stack.pop_back();

            auto relaxation = simplex_maximize(node.a, node.b, relaxation_objective, eps);
            if (relaxation.is_infeasible()) continue;

            if (relaxation.is_unbounded()) {
                std::vector<T> feasible_variables;
                if (has_positive_direction(node) && find_integer_feasible(node, feasible_variables)) {
                    return make_unbounded_result(feasible_variables);
                }
                continue;
            }

            if (can_prune_by_bound(relaxation.objective_value)) continue;

            std::vector<T> integer_variables;
            if (round_solution(relaxation.variables, integer_variables)) {
                update_incumbent(integer_variables);
                continue;
            }

            int variable = find_fractional_variable(relaxation.variables);
            if (variable == -1) continue;
            push_branches(stack, node, variable, relaxation.variables[variable]);
        }

        if (!has_incumbent) return make_infeasible_result();
        return make_optimal_result();
    }
};

}  // namespace detail

template <class T>
IntegerLpResult<T> integer_lp_maximize(const std::vector<std::vector<T>>& a,
                                       const std::vector<T>& b, const std::vector<T>& c,
                                       long double eps = 1e-10L) {
    static_assert(std::is_integral_v<T> && std::is_signed_v<T>,
                  "integer_lp requires a signed integer type");
    assert(int(a.size()) == int(b.size()));
    for (const auto& row : a) assert(int(row.size()) == int(c.size()));
    assert(eps > 0);

    detail::IntegerLpSolver<T> solver(a, b, c, true, eps);
    return solver.solve();
}

template <class T>
IntegerLpResult<T> integer_lp_minimize(const std::vector<std::vector<T>>& a,
                                       const std::vector<T>& b, const std::vector<T>& c,
                                       long double eps = 1e-10L) {
    static_assert(std::is_integral_v<T> && std::is_signed_v<T>,
                  "integer_lp requires a signed integer type");
    assert(int(a.size()) == int(b.size()));
    for (const auto& row : a) assert(int(row.size()) == int(c.size()));
    assert(eps > 0);

    detail::IntegerLpSolver<T> solver(a, b, c, false, eps);
    return solver.solve();
}

template <class T>
IntegerLpResult<T> integer_lp(const std::vector<std::vector<T>>& a, const std::vector<T>& b,
                              const std::vector<T>& c, long double eps = 1e-10L) {
    return integer_lp_maximize(a, b, c, eps);
}

}  // namespace opt
}  // namespace m1une

#endif  // M1UNE_OPTIMIZATION_INTEGER_LP_HPP

#line 1 "optimization/integer_lp.hpp"



#include <algorithm>
#include <cassert>
#include <cmath>
#include <limits>
#include <type_traits>
#include <vector>

#line 1 "optimization/simplex.hpp"



#line 7 "optimization/simplex.hpp"
#include <utility>
#line 9 "optimization/simplex.hpp"

namespace m1une {
namespace opt {

enum class SimplexStatus {
    Optimal,
    Infeasible,
    Unbounded,
};

template <class T>
struct SimplexResult {
    SimplexStatus status;
    T objective_value;
    std::vector<T> variables;

    bool is_optimal() const { return status == SimplexStatus::Optimal; }
    bool is_infeasible() const { return status == SimplexStatus::Infeasible; }
    bool is_unbounded() const { return status == SimplexStatus::Unbounded; }
};

namespace detail {

template <class T>
T simplex_abs(T x) {
    return x < T() ? -x : x;
}

template <class T>
struct SimplexTableau {
    int constraint_count;
    int variable_count;
    T eps;
    std::vector<int> basis;
    std::vector<int> nonbasis;
    std::vector<std::vector<T>> table;

    SimplexTableau(const std::vector<std::vector<T>>& a, const std::vector<T>& b,
                   const std::vector<T>& c, T epsilon)
        : constraint_count(int(b.size())),
          variable_count(int(c.size())),
          eps(epsilon),
          basis(constraint_count),
          nonbasis(variable_count + 1),
          table(constraint_count + 2, std::vector<T>(variable_count + 2, T())) {
        for (int i = 0; i < constraint_count; i++) {
            for (int j = 0; j < variable_count; j++) table[i][j] = a[i][j];
        }
        for (int i = 0; i < constraint_count; i++) {
            basis[i] = variable_count + i;
            table[i][artificial_col()] = T(-1);
            table[i][rhs_col()] = b[i];
        }
        for (int j = 0; j < variable_count; j++) {
            nonbasis[j] = j;
            table[objective_row()][j] = -c[j];
        }
        nonbasis[artificial_col()] = artificial_id();
        table[auxiliary_row()][artificial_col()] = T(1);
    }

    int objective_row() const { return constraint_count; }
    int auxiliary_row() const { return constraint_count + 1; }
    int artificial_col() const { return variable_count; }
    int rhs_col() const { return variable_count + 1; }
    int artificial_id() const { return -1; }

    T normalize(T x) const {
        return simplex_abs(x) <= eps ? T() : x;
    }

    bool less_with_tie(int row, int lhs, int rhs) const {
        if (table[row][lhs] < table[row][rhs] - eps) return true;
        if (table[row][rhs] < table[row][lhs] - eps) return false;
        return nonbasis[lhs] < nonbasis[rhs];
    }

    bool better_leaving_row(int lhs, int rhs, int entering_col) const {
        T lhs_ratio = table[lhs][rhs_col()] / table[lhs][entering_col];
        T rhs_ratio = table[rhs][rhs_col()] / table[rhs][entering_col];
        if (lhs_ratio < rhs_ratio - eps) return true;
        if (rhs_ratio < lhs_ratio - eps) return false;
        return basis[lhs] < basis[rhs];
    }

    void pivot(int leaving_row, int entering_col) {
        T inverse = T(1) / table[leaving_row][entering_col];
        for (int i = 0; i < constraint_count + 2; i++) {
            if (i == leaving_row) continue;
            for (int j = 0; j < variable_count + 2; j++) {
                if (j == entering_col) continue;
                table[i][j] -= table[leaving_row][j] * table[i][entering_col] * inverse;
            }
        }
        for (int j = 0; j < variable_count + 2; j++) {
            if (j != entering_col) table[leaving_row][j] *= inverse;
        }
        for (int i = 0; i < constraint_count + 2; i++) {
            if (i != leaving_row) table[i][entering_col] *= -inverse;
        }
        table[leaving_row][entering_col] = inverse;
        std::swap(basis[leaving_row], nonbasis[entering_col]);
    }

    bool run_simplex(int row) {
        while (true) {
            int entering_col = -1;
            for (int j = 0; j <= variable_count; j++) {
                if (nonbasis[j] == artificial_id()) continue;
                if (entering_col == -1 || less_with_tie(row, j, entering_col)) entering_col = j;
            }
            if (entering_col == -1 || table[row][entering_col] >= -eps) return true;

            int leaving_row = -1;
            for (int i = 0; i < constraint_count; i++) {
                if (table[i][entering_col] <= eps) continue;
                if (leaving_row == -1 || better_leaving_row(i, leaving_row, entering_col)) {
                    leaving_row = i;
                }
            }
            if (leaving_row == -1) return false;
            pivot(leaving_row, entering_col);
        }
    }

    bool make_feasible() {
        int leaving_row = 0;
        for (int i = 1; i < constraint_count; i++) {
            if (table[i][rhs_col()] < table[leaving_row][rhs_col()]) leaving_row = i;
        }
        if (constraint_count == 0 || table[leaving_row][rhs_col()] >= -eps) return true;

        pivot(leaving_row, artificial_col());
        if (!run_simplex(auxiliary_row())) return false;
        if (table[auxiliary_row()][rhs_col()] < -eps) return false;

        for (int i = 0; i < constraint_count; i++) {
            if (basis[i] != artificial_id()) continue;
            int entering_col = -1;
            for (int j = 0; j <= variable_count; j++) {
                if (nonbasis[j] == artificial_id()) continue;
                if (simplex_abs(table[i][j]) <= eps) continue;
                if (entering_col == -1 || nonbasis[j] < nonbasis[entering_col]) entering_col = j;
            }
            if (entering_col != -1) pivot(i, entering_col);
        }
        return true;
    }

    SimplexStatus solve(std::vector<T>& variables, T& objective_value) {
        if (!make_feasible()) return SimplexStatus::Infeasible;
        if (!run_simplex(objective_row())) return SimplexStatus::Unbounded;

        variables.assign(variable_count, T());
        for (int i = 0; i < constraint_count; i++) {
            if (0 <= basis[i] && basis[i] < variable_count) {
                variables[basis[i]] = normalize(table[i][rhs_col()]);
            }
        }
        objective_value = normalize(table[objective_row()][rhs_col()]);
        return SimplexStatus::Optimal;
    }
};

}  // namespace detail

template <class T>
SimplexResult<T> simplex_maximize(const std::vector<std::vector<T>>& a, const std::vector<T>& b,
                                  const std::vector<T>& c, T eps = T(1e-10)) {
    static_assert(std::is_floating_point_v<T>, "simplex requires a floating-point type");
    assert(int(a.size()) == int(b.size()));
    for (const auto& row : a) assert(int(row.size()) == int(c.size()));
    assert(eps > T());

    SimplexResult<T> result;
    result.status = SimplexStatus::Infeasible;
    result.objective_value = std::numeric_limits<T>::quiet_NaN();
    result.variables.assign(c.size(), T());

    detail::SimplexTableau<T> solver(a, b, c, eps);
    result.status = solver.solve(result.variables, result.objective_value);
    if (result.status == SimplexStatus::Infeasible) {
        result.objective_value = std::numeric_limits<T>::quiet_NaN();
    } else if (result.status == SimplexStatus::Unbounded) {
        result.objective_value = std::numeric_limits<T>::infinity();
    }
    return result;
}

template <class T>
SimplexResult<T> simplex_minimize(const std::vector<std::vector<T>>& a, const std::vector<T>& b,
                                  const std::vector<T>& c, T eps = T(1e-10)) {
    std::vector<T> negated = c;
    for (T& x : negated) x = -x;
    auto result = simplex_maximize(a, b, negated, eps);
    if (result.status == SimplexStatus::Optimal) {
        result.objective_value = -result.objective_value;
    } else if (result.status == SimplexStatus::Unbounded) {
        result.objective_value = -std::numeric_limits<T>::infinity();
    }
    return result;
}

template <class T>
SimplexResult<T> simplex(const std::vector<std::vector<T>>& a, const std::vector<T>& b,
                         const std::vector<T>& c, T eps = T(1e-10)) {
    return simplex_maximize(a, b, c, eps);
}

}  // namespace opt
}  // namespace m1une


#line 12 "optimization/integer_lp.hpp"

namespace m1une {
namespace opt {

enum class IntegerLpStatus {
    Optimal,
    Infeasible,
    Unbounded,
};

template <class T>
struct IntegerLpResult {
    IntegerLpStatus status;
    T objective_value;
    std::vector<T> variables;

    bool is_optimal() const { return status == IntegerLpStatus::Optimal; }
    bool is_infeasible() const { return status == IntegerLpStatus::Infeasible; }
    bool is_unbounded() const { return status == IntegerLpStatus::Unbounded; }
};

namespace detail {

template <class T>
struct IntegerLpSolver {
    using Real = long double;

    struct Node {
        std::vector<std::vector<Real>> a;
        std::vector<Real> b;
    };

    int variable_count;
    bool maximize;
    Real eps;
    std::vector<T> objective;
    std::vector<Real> relaxation_objective;
    Node initial_node;

    bool has_incumbent = false;
    T best_value = T();
    std::vector<T> best_variables;

    IntegerLpSolver(const std::vector<std::vector<T>>& a, const std::vector<T>& b,
                    const std::vector<T>& c, bool is_maximize, Real epsilon)
        : variable_count(int(c.size())),
          maximize(is_maximize),
          eps(epsilon),
          objective(c),
          relaxation_objective(c.size(), Real()),
          initial_node() {
        initial_node.a.assign(a.size(), std::vector<Real>(variable_count, Real()));
        initial_node.b.assign(b.size(), Real());
        for (int i = 0; i < int(a.size()); i++) {
            for (int j = 0; j < variable_count; j++) initial_node.a[i][j] = Real(a[i][j]);
            initial_node.b[i] = Real(b[i]);
        }
        Real sign = maximize ? Real(1) : Real(-1);
        for (int j = 0; j < variable_count; j++) relaxation_objective[j] = sign * Real(c[j]);
    }

    Real abs_value(Real x) const {
        return x < Real() ? -x : x;
    }

    bool better_value(T lhs, T rhs) const {
        return maximize ? lhs > rhs : lhs < rhs;
    }

    bool can_prune_by_bound(Real relaxation_value) const {
        if (!has_incumbent) return false;
        Real signed_best = maximize ? Real(best_value) : -Real(best_value);
        return relaxation_value <= signed_best + eps;
    }

    T evaluate(const std::vector<T>& variables) const {
        T result = T();
        for (int i = 0; i < variable_count; i++) result += objective[i] * variables[i];
        return result;
    }

    bool round_solution(const std::vector<Real>& real_variables, std::vector<T>& variables) const {
        variables.assign(variable_count, T());
        for (int i = 0; i < variable_count; i++) {
            Real value = real_variables[i];
            if (value < -eps) return false;
            Real rounded = std::round(value);
            if (abs_value(value - rounded) > eps) return false;
            variables[i] = static_cast<T>(rounded);
        }
        return true;
    }

    int find_fractional_variable(const std::vector<Real>& real_variables) const {
        int result = -1;
        Real best_distance = eps;
        for (int i = 0; i < variable_count; i++) {
            Real value = real_variables[i];
            Real rounded = std::round(value);
            Real distance = abs_value(value - rounded);
            if (distance > best_distance) {
                best_distance = distance;
                result = i;
            }
        }
        return result;
    }

    Node with_upper_bound(const Node& node, int variable, T bound) const {
        Node result = node;
        result.a.emplace_back(variable_count, Real());
        result.a.back()[variable] = Real(1);
        result.b.push_back(Real(bound));
        return result;
    }

    Node with_lower_bound(const Node& node, int variable, T bound) const {
        Node result = node;
        result.a.emplace_back(variable_count, Real());
        result.a.back()[variable] = Real(-1);
        result.b.push_back(-Real(bound));
        return result;
    }

    void push_branches(std::vector<Node>& stack, const Node& node, int variable, Real value) const {
        Real floor_value = std::floor(value);
        Real ceil_value = std::ceil(value);
        T upper_bound = static_cast<T>(floor_value);
        T lower_bound = static_cast<T>(ceil_value);

        bool has_upper_branch = upper_bound >= T();
        bool prefer_lower_branch = relaxation_objective[variable] >= -eps;

        if (prefer_lower_branch) {
            if (has_upper_branch) stack.push_back(with_upper_bound(node, variable, upper_bound));
            stack.push_back(with_lower_bound(node, variable, lower_bound));
        } else {
            stack.push_back(with_lower_bound(node, variable, lower_bound));
            if (has_upper_branch) stack.push_back(with_upper_bound(node, variable, upper_bound));
        }
    }

    bool has_positive_direction(const Node& node) const {
        std::vector<std::vector<Real>> direction_a = node.a;
        std::vector<Real> direction_b(node.b.size(), Real());

        std::vector<Real> objective_row(variable_count, Real());
        for (int i = 0; i < variable_count; i++) objective_row[i] = -relaxation_objective[i];
        direction_a.push_back(objective_row);
        direction_b.push_back(Real(-1));

        std::vector<Real> zero_objective(variable_count, Real());
        auto result = simplex_maximize(direction_a, direction_b, zero_objective, eps);
        return result.is_optimal();
    }

    bool find_integer_feasible(const Node& start, std::vector<T>& feasible_variables) const {
        std::vector<Node> stack;
        stack.push_back(start);
        std::vector<Real> zero_objective(variable_count, Real());

        while (!stack.empty()) {
            Node node = stack.back();
            stack.pop_back();

            auto relaxation = simplex_maximize(node.a, node.b, zero_objective, eps);
            if (relaxation.is_infeasible()) continue;
            if (relaxation.is_unbounded()) continue;

            if (round_solution(relaxation.variables, feasible_variables)) return true;

            int variable = find_fractional_variable(relaxation.variables);
            if (variable == -1) continue;
            push_branches(stack, node, variable, relaxation.variables[variable]);
        }
        return false;
    }

    void update_incumbent(const std::vector<T>& variables) {
        T value = evaluate(variables);
        if (!has_incumbent || better_value(value, best_value)) {
            has_incumbent = true;
            best_value = value;
            best_variables = variables;
        }
    }

    IntegerLpResult<T> make_infeasible_result() const {
        IntegerLpResult<T> result;
        result.status = IntegerLpStatus::Infeasible;
        result.objective_value = T();
        result.variables.assign(variable_count, T());
        return result;
    }

    IntegerLpResult<T> make_unbounded_result(const std::vector<T>& variables) const {
        IntegerLpResult<T> result;
        result.status = IntegerLpStatus::Unbounded;
        result.objective_value =
            maximize ? std::numeric_limits<T>::max() : std::numeric_limits<T>::lowest();
        result.variables = variables;
        return result;
    }

    IntegerLpResult<T> make_optimal_result() const {
        IntegerLpResult<T> result;
        result.status = IntegerLpStatus::Optimal;
        result.objective_value = best_value;
        result.variables = best_variables;
        return result;
    }

    IntegerLpResult<T> solve() {
        std::vector<Node> stack;
        stack.push_back(initial_node);

        while (!stack.empty()) {
            Node node = stack.back();
            stack.pop_back();

            auto relaxation = simplex_maximize(node.a, node.b, relaxation_objective, eps);
            if (relaxation.is_infeasible()) continue;

            if (relaxation.is_unbounded()) {
                std::vector<T> feasible_variables;
                if (has_positive_direction(node) && find_integer_feasible(node, feasible_variables)) {
                    return make_unbounded_result(feasible_variables);
                }
                continue;
            }

            if (can_prune_by_bound(relaxation.objective_value)) continue;

            std::vector<T> integer_variables;
            if (round_solution(relaxation.variables, integer_variables)) {
                update_incumbent(integer_variables);
                continue;
            }

            int variable = find_fractional_variable(relaxation.variables);
            if (variable == -1) continue;
            push_branches(stack, node, variable, relaxation.variables[variable]);
        }

        if (!has_incumbent) return make_infeasible_result();
        return make_optimal_result();
    }
};

}  // namespace detail

template <class T>
IntegerLpResult<T> integer_lp_maximize(const std::vector<std::vector<T>>& a,
                                       const std::vector<T>& b, const std::vector<T>& c,
                                       long double eps = 1e-10L) {
    static_assert(std::is_integral_v<T> && std::is_signed_v<T>,
                  "integer_lp requires a signed integer type");
    assert(int(a.size()) == int(b.size()));
    for (const auto& row : a) assert(int(row.size()) == int(c.size()));
    assert(eps > 0);

    detail::IntegerLpSolver<T> solver(a, b, c, true, eps);
    return solver.solve();
}

template <class T>
IntegerLpResult<T> integer_lp_minimize(const std::vector<std::vector<T>>& a,
                                       const std::vector<T>& b, const std::vector<T>& c,
                                       long double eps = 1e-10L) {
    static_assert(std::is_integral_v<T> && std::is_signed_v<T>,
                  "integer_lp requires a signed integer type");
    assert(int(a.size()) == int(b.size()));
    for (const auto& row : a) assert(int(row.size()) == int(c.size()));
    assert(eps > 0);

    detail::IntegerLpSolver<T> solver(a, b, c, false, eps);
    return solver.solve();
}

template <class T>
IntegerLpResult<T> integer_lp(const std::vector<std::vector<T>>& a, const std::vector<T>& b,
                              const std::vector<T>& c, long double eps = 1e-10L) {
    return integer_lp_maximize(a, b, c, eps);
}

}  // namespace opt
}  // namespace m1une

m1une's library

This documentation is automatically generated by online-judge-tools/verification-helper

Integer Linear Programming
(optimization/integer_lp.hpp)

Overview

Interface

Functions

Example

Depends on

Required by

Verified with

Code

Integer Linear Programming (optimization/integer_lp.hpp)

Overview

Interface

Functions

Example

Depends on

Required by

Verified with

Code

Integer Linear Programming
(optimization/integer_lp.hpp)