m1une's library

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

View on GitHub

:heavy_check_mark: Project Selection
(optimization/project_selection.hpp)

Overview

ProjectSelection<T> chooses a subset of binary projects that maximizes total gain. Besides a gain or cost for each individual choice, it supports implication penalties, hard implications, forced choices, and several common group rewards.

The problem is reduced to a minimum s-t cut and solved with Dinic’s algorithm. This is substantially faster and more predictable than a general integer programming solver for objectives that fit the supported forms.

For ordered variables with more than two choices and supermodular pairwise gain tables, use KProjectSelection<T>.

Project i is selected when result.selected[i] is true. Gains increase the objective and penalties decrease it. A cost of c for selecting project i can therefore be written as add_gain(i, -c).

Supported Objective Terms

All gains and penalties are additive.

Call Contribution to the objective
add_gain(i, selected_gain) Adds selected_gain if i is selected.
add_gain(i, selected_gain, unselected_gain) Adds the gain corresponding to the choice for i.
add_penalty(i, j, penalty) Subtracts penalty exactly when i is selected and j is unselected.
add_penalty_if_different(i, j, penalty) Subtracts penalty when exactly one of i and j is selected.
add_gain_if_same(i, j, gain) Adds gain when i and j make the same choice.
add_gain_if_all_selected(projects, gain) Adds gain when every listed project is selected.
add_gain_if_all_unselected(projects, gain) Adds gain when every listed project is unselected.

Domain of Gains and Penalties

The unary gains passed to either overload of add_gain may be negative, zero, or positive. In particular, a negative gain represents a cost for that choice. Both selected_gain and unselected_gain may be negative, and the optimal max_gain may also be negative when every feasible selection has a negative total gain.

The other numeric arguments have a restricted domain because they become flow capacities:

Argument Required domain
add_gain(i, selected_gain) selected_gain may be any value representable by T.
add_gain(i, selected_gain, unselected_gain) Both gains may be any values representable by T.
add_penalty, add_penalty_if_different penalty >= 0.
add_gain_if_same gain >= 0.
add_gain_if_all_selected, add_gain_if_all_unselected gain >= 0.

The group methods treat an empty list as vacuously satisfying the condition, so their non-negative gain is always added. All accumulated values and intermediate differences must additionally satisfy the range requirements in the Numeric Requirements section.

add_penalty(i, j, penalty) is the low-level finite implication primitive. It is useful for statements such as “choosing i without choosing j loses 20.” It does not forbid that combination; use add_hard_implication when it must be impossible.

Not every Boolean objective can be represented by one minimum cut. In particular, a positive reward for “at least one is selected” or a penalty for “both are selected” is not generally supported without reformulating the problem. The methods above are exactly the forms provided by this interface.

Reduction to Minimum Cut

The source side of the cut represents selected projects, and the sink side represents unselected projects.

Thus every cut has cost equal to a constant minus the modeled gain. A minimum cut therefore gives a maximum-gain selection. If several selections have the same maximum gain, solve() may return any one of them.

Hard Constraints

Call Constraint
add_hard_implication(i, j) If i is selected, j must also be selected.
force_selected(i) Project i must be selected.
force_unselected(i) Project i must be unselected.

Hard constraints use a capacity larger than the sum of every finite capacity. This value is computed only inside solve(), so constraints remain truly hard even if more gains or penalties are added later.

Contradictory hard constraints make the result infeasible. Check result.is_feasible() before reading max_gain or interpreting selected.

Result

ProjectSelectionResult<T> contains:

Member / Method Type / Signature Meaning
feasible bool Whether all hard constraints can be satisfied.
max_gain T Maximum total gain. It is meaningful only when feasible.
selected std::vector<bool> One optimal choice for the original projects. Auxiliary vertices are not included.
is_feasible bool is_feasible() const Returns feasible.

Calling solve() does not mutate the model, so it may be called repeatedly.

Methods and Complexity

Let N be the number of original projects plus internally created auxiliary vertices, and let M be the number of generated flow edges.

Method Complexity
Constructor, size $O(1)$
add_gain, add_penalty, add_penalty_if_different, add_gain_if_same, hard single-project methods Amortized $O(1)$
add_gain_if_all_selected, add_gain_if_all_unselected $O(K)$ for a group of size K
solve General-case $O(N^2 M)$ time and $O(N + M)$ memory

A group-reward call creates one auxiliary vertex when its group has at least three entries. solve() uses the repository’s MaxFlow<T> implementation.

Numeric Requirements

T must be a signed integral type; long long is recommended. Every intermediate gain difference, the sum of finite capacities, the hard capacity, and the final answer must fit in T. These range requirements are checked with assertions where practical.

Example

Suppose each selected project earns its listed gain, project 0 requires project 1, and completing projects 1 and 2 together earns a bonus.

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

int main() {
    m1une::opt::ProjectSelection<long long> solver(3);
    solver.add_gain(0, 10);
    solver.add_gain(1, -3);  // Selecting project 1 costs 3.
    solver.add_gain(2, 4);
    solver.add_hard_implication(0, 1);
    solver.add_gain_if_all_selected(std::vector<int>{1, 2}, 5);

    auto result = solver.solve();
    if (!result.is_feasible()) return 0;

    std::cout << result.max_gain << "\n";  // 16
    for (int i = 0; i < solver.size(); i++) {
        if (result.selected[i]) std::cout << i << "\n";
    }
}

Depends on

Required by

Verified with

Code

#ifndef M1UNE_OPTIMIZATION_PROJECT_SELECTION_HPP
#define M1UNE_OPTIMIZATION_PROJECT_SELECTION_HPP 1

#include <cassert>
#include <limits>
#include <type_traits>
#include <utility>
#include <vector>

#include "../graph/flow/max_flow.hpp"

namespace m1une {
namespace opt {

template <class T>
struct ProjectSelectionResult {
    bool feasible;
    T max_gain;
    std::vector<bool> selected;

    bool is_feasible() const {
        return feasible;
    }
};

template <class T>
class ProjectSelection {
    static_assert(std::is_integral_v<T> && std::is_signed_v<T>);

    struct Arc {
        int from;
        int to;
        T cap;
    };

    static constexpr int source = -1;
    static constexpr int sink = -2;

    int _project_count;
    int _vertex_count;
    T _offset = T();
    T _finite_cap_sum = T();
    std::vector<Arc> _arcs;
    std::vector<std::pair<int, int>> _hard_arcs;

    void assert_project(int project) const {
        (void)project;
        assert(0 <= project && project < _project_count);
    }

    void assert_vertex(int vertex) const {
        (void)vertex;
        assert(0 <= vertex && vertex < _vertex_count);
    }

    void add_offset(T value) {
        if (value > T()) {
            assert(_offset <= std::numeric_limits<T>::max() - value);
        } else if (value < T()) {
            assert(_offset >= std::numeric_limits<T>::lowest() - value);
        }
        _offset += value;
    }

    T nonnegative_difference(T large, T small) const {
        assert(small <= large);
        if (small < T()) {
            assert(large <= std::numeric_limits<T>::max() + small);
        }
        return large - small;
    }

    void add_arc(int from, int to, T cap) {
        assert(cap >= T());
        if (from == to) return;
        assert(cap <= std::numeric_limits<T>::max() - _finite_cap_sum);
        _finite_cap_sum += cap;
        _arcs.push_back(Arc{from, to, cap});
    }

    void add_hard_arc(int from, int to) {
        if (from == to) return;
        _hard_arcs.emplace_back(from, to);
    }

    void add_vertex_gain(int vertex, T gain_if_selected, T gain_if_unselected) {
        assert_vertex(vertex);
        if (gain_if_selected >= gain_if_unselected) {
            add_offset(gain_if_selected);
            add_arc(source, vertex,
                    nonnegative_difference(gain_if_selected, gain_if_unselected));
        } else {
            add_offset(gain_if_unselected);
            add_arc(vertex, sink,
                    nonnegative_difference(gain_if_unselected, gain_if_selected));
        }
    }

    int add_auxiliary_vertex() {
        return _vertex_count++;
    }

   public:
    ProjectSelection() : ProjectSelection(0) {}

    explicit ProjectSelection(int project_count)
        : _project_count(project_count), _vertex_count(project_count) {
        assert(project_count >= 0);
    }

    int size() const {
        return _project_count;
    }

    void add_gain(int project, T gain_if_selected) {
        add_gain(project, gain_if_selected, T());
    }

    void add_gain(int project, T gain_if_selected, T gain_if_unselected) {
        assert_project(project);
        add_vertex_gain(project, gain_if_selected, gain_if_unselected);
    }

    void add_penalty(int selected_project, int unselected_project, T penalty) {
        assert_project(selected_project);
        assert_project(unselected_project);
        add_arc(selected_project, unselected_project, penalty);
    }

    void add_penalty_if_different(int project_a, int project_b, T penalty) {
        assert_project(project_a);
        assert_project(project_b);
        add_arc(project_a, project_b, penalty);
        add_arc(project_b, project_a, penalty);
    }

    void add_gain_if_same(int project_a, int project_b, T gain) {
        assert(gain >= T());
        add_offset(gain);
        add_penalty_if_different(project_a, project_b, gain);
    }

    void add_hard_implication(int selected_project, int required_project) {
        assert_project(selected_project);
        assert_project(required_project);
        add_hard_arc(selected_project, required_project);
    }

    void force_selected(int project) {
        assert_project(project);
        add_hard_arc(source, project);
    }

    void force_unselected(int project) {
        assert_project(project);
        add_hard_arc(project, sink);
    }

    void add_gain_if_all_selected(const std::vector<int>& projects, T gain) {
        assert(gain >= T());
        for (int project : projects) assert_project(project);
        if (projects.empty()) {
            add_offset(gain);
            return;
        }
        if (projects.size() == 1) {
            add_vertex_gain(projects[0], gain, T());
            return;
        }
        if (projects.size() == 2) {
            add_vertex_gain(projects[0], gain, T());
            add_arc(projects[0], projects[1], gain);
            return;
        }

        int auxiliary = add_auxiliary_vertex();
        add_vertex_gain(auxiliary, gain, T());
        for (int project : projects) add_hard_arc(auxiliary, project);
    }

    void add_gain_if_all_unselected(const std::vector<int>& projects, T gain) {
        assert(gain >= T());
        for (int project : projects) assert_project(project);
        if (projects.empty()) {
            add_offset(gain);
            return;
        }
        if (projects.size() == 1) {
            add_vertex_gain(projects[0], T(), gain);
            return;
        }
        if (projects.size() == 2) {
            add_vertex_gain(projects[0], T(), gain);
            add_arc(projects[1], projects[0], gain);
            return;
        }

        int auxiliary = add_auxiliary_vertex();
        add_vertex_gain(auxiliary, T(), gain);
        for (int project : projects) add_hard_arc(project, auxiliary);
    }

    ProjectSelectionResult<T> solve() const {
        int s = _vertex_count;
        int t = s + 1;
        flow::MaxFlow<T> max_flow(_vertex_count + 2);

        auto vertex_id = [&](int vertex) {
            if (vertex == source) return s;
            if (vertex == sink) return t;
            return vertex;
        };

        for (const auto& arc : _arcs) {
            max_flow.add_edge(vertex_id(arc.from), vertex_id(arc.to), arc.cap);
        }

        T hard_cap = T();
        if (!_hard_arcs.empty()) {
            assert(_finite_cap_sum < std::numeric_limits<T>::max());
            hard_cap = _finite_cap_sum + T(1);
            for (auto [from, to] : _hard_arcs) {
                max_flow.add_edge(vertex_id(from), vertex_id(to), hard_cap);
            }
        }

        T cut_cost =
            _hard_arcs.empty() ? max_flow.max_flow(s, t) : max_flow.max_flow(s, t, hard_cap);
        ProjectSelectionResult<T> result;
        result.feasible = _hard_arcs.empty() || cut_cost < hard_cap;
        result.max_gain = T();
        result.selected.assign(_project_count, false);
        if (!result.feasible) return result;

        assert(_offset >= std::numeric_limits<T>::lowest() + cut_cost);
        result.max_gain = _offset - cut_cost;
        auto source_side = max_flow.min_cut(s);
        for (int project = 0; project < _project_count; project++) {
            result.selected[project] = source_side[project];
        }
        return result;
    }
};

}  // namespace opt
}  // namespace m1une

#endif  // M1UNE_OPTIMIZATION_PROJECT_SELECTION_HPP
#line 1 "optimization/project_selection.hpp"



#include <cassert>
#include <limits>
#include <type_traits>
#include <utility>
#include <vector>

#line 1 "graph/flow/max_flow.hpp"



#include <algorithm>
#line 6 "graph/flow/max_flow.hpp"
#include <cstddef>
#line 9 "graph/flow/max_flow.hpp"

namespace m1une {
namespace flow {

template <class Cap>
struct MaxFlow {
    struct Edge {
        int from;
        int to;
        Cap cap;
        Cap flow;
    };

   private:
    struct InternalEdge {
        int to;
        int rev;
        Cap cap;
    };

    struct Position {
        int from;
        int edge;
    };

    int _n;
    std::vector<Position> _pos;
    std::vector<std::vector<InternalEdge>> _g;

    Cap highest_label_preflow_push(int s, int t) {
        const int dead = 2 * _n;
        const int unreachable = _n + 1;
        std::vector<Cap> excess(_n, Cap(0));
        std::vector<int> state(8 * std::size_t(_n) + 2);
        int* height = state.data();
        int* height_count = height + _n;
        int* current = height_count + dead + 1;
        int* queue = current + _n;
        int* next = queue + _n;
        int* bucket_head = next + _n;
        std::vector<char> active(_n, false);
        int highest = -1;
        long long work = 0;
        const long long arc_count =
            2LL * static_cast<long long>(_pos.size());
        const long long work_limit = std::max(1LL, 4 * arc_count + _n);

        auto activate = [&](int v) {
            if (v == s || v == t || active[v] || excess[v] == Cap(0) ||
                height[v] >= dead) {
                return;
            }
            active[v] = true;
            next[v] = bucket_head[height[v]];
            bucket_head[height[v]] = v;
            highest = std::max(highest, height[v]);
        };

        auto rebuild_buckets = [&]() {
            std::fill(bucket_head, bucket_head + dead + 1, -1);
            std::fill(active.begin(), active.end(), false);
            highest = -1;
            for (int v = 0; v < _n; v++) activate(v);
        };

        auto global_relabel = [&]() {
            std::fill(height, height + _n, unreachable);
            std::fill(height_count, height_count + dead + 1, 0);
            std::fill(current, current + _n, 0);
            int head = 0;
            int tail = 0;
            height[t] = 0;
            height[s] = _n;
            queue[tail++] = t;
            while (head != tail) {
                int v = queue[head++];
                for (const auto& e : _g[v]) {
                    if (e.to == s || height[e.to] != unreachable) continue;
                    const auto& reverse = _g[e.to][e.rev];
                    if (reverse.cap == Cap(0)) continue;
                    height[e.to] = height[v] + 1;
                    queue[tail++] = e.to;
                }
            }
            for (int v = 0; v < _n; v++) height_count[height[v]]++;
            rebuild_buckets();
            work = 0;
        };

        auto gap = [&](int empty_height) {
            for (int v = 0; v < _n; v++) {
                if (v == s || v == t || height[v] <= empty_height ||
                    height[v] >= _n) {
                    continue;
                }
                height_count[height[v]]--;
                height[v] = unreachable;
                height_count[height[v]]++;
                current[v] = 0;
            }
            rebuild_buckets();
        };

        auto relabel = [&](int v) -> bool {
            int old_height = height[v];
            int new_height = dead;
            work += int(_g[v].size());
            for (const auto& e : _g[v]) {
                if (e.cap != Cap(0)) {
                    new_height = std::min(new_height, height[e.to] + 1);
                }
            }
            height_count[old_height]--;
            height[v] = std::min(new_height, dead);
            height_count[height[v]]++;
            current[v] = 0;
            if (old_height < _n && height_count[old_height] == 0) {
                gap(old_height);
                return true;
            }
            return false;
        };

        auto push = [&](int v, InternalEdge& e) {
            Cap sent = std::min(excess[v], e.cap);
            bool was_zero = excess[e.to] == Cap(0);
            e.cap -= sent;
            _g[e.to][e.rev].cap += sent;
            excess[v] -= sent;
            excess[e.to] += sent;
            if (was_zero) activate(e.to);
        };

        auto discharge = [&](int v) {
            while (excess[v] != Cap(0) && height[v] < dead) {
                if (current[v] == int(_g[v].size())) {
                    if (relabel(v)) return;
                    continue;
                }
                auto& e = _g[v][current[v]];
                work++;
                if (e.cap != Cap(0) && height[v] == height[e.to] + 1) {
                    push(v, e);
                } else {
                    current[v]++;
                }
            }
            activate(v);
        };

        for (auto& e : _g[s]) {
            if (e.to == s || e.cap == Cap(0)) continue;
            Cap sent = e.cap;
            e.cap = Cap(0);
            _g[e.to][e.rev].cap += sent;
            excess[e.to] += sent;
        }
        global_relabel();

        while (highest >= 0) {
            if (bucket_head[highest] == -1) {
                highest--;
                continue;
            }
            int v = bucket_head[highest];
            bucket_head[highest] = next[v];
            if (!active[v] || height[v] != highest) continue;
            active[v] = false;
            discharge(v);
            if (work >= work_limit) global_relabel();
        }
        return excess[t];
    }

   public:
    MaxFlow() : MaxFlow(0) {}

    explicit MaxFlow(int n) : _n(n), _g(n) {
        assert(0 <= n);
    }

    int size() const {
        return _n;
    }

    int edge_count() const {
        return int(_pos.size());
    }

    void reserve_edges(int edge_count) {
        assert(0 <= edge_count);
        _pos.reserve(edge_count);
        if (_n == 0 || edge_count == 0 ||
            2 * std::size_t(edge_count) < std::size_t(_n)) {
            return;
        }
        const std::size_t average_degree =
            (3 * std::size_t(edge_count) + std::size_t(_n) - 1)
            / std::size_t(_n);
        for (auto& edges : _g) edges.reserve(average_degree);
    }

    void reserve_edges(int edge_count, const std::vector<int>& degrees) {
        assert(0 <= edge_count);
        assert(int(degrees.size()) == _n);
        _pos.reserve(edge_count);
        for (int v = 0; v < _n; v++) {
            assert(0 <= degrees[v]);
            _g[v].reserve(degrees[v]);
        }
    }

    int add_edge(int from, int to, Cap cap) {
        assert(0 <= from && from < _n);
        assert(0 <= to && to < _n);
        assert(Cap(0) <= cap);
        int id = int(_pos.size());
        int from_id = int(_g[from].size());
        int to_id = int(_g[to].size());
        if (from == to) to_id++;
        _pos.push_back(Position{from, from_id});
        _g[from].push_back(InternalEdge{to, to_id, cap});
        _g[to].push_back(InternalEdge{from, from_id, Cap(0)});
        return id;
    }

    int add_undirected_edge(int first, int second, Cap cap) {
        static_assert(std::numeric_limits<Cap>::is_signed);
        assert(0 <= first && first < _n);
        assert(0 <= second && second < _n);
        assert(Cap(0) <= cap);
        assert(cap <= std::numeric_limits<Cap>::max() / Cap(2));
        int id = int(_pos.size());
        int first_id = int(_g[first].size());
        int second_id = int(_g[second].size());
        if (first == second) second_id++;
        _pos.push_back(Position{first, ~first_id});
        _g[first].push_back(InternalEdge{second, second_id, cap});
        _g[second].push_back(InternalEdge{first, first_id, cap});
        return id;
    }

    Edge get_edge(int i) const {
        assert(0 <= i && i < int(_pos.size()));
        const auto& position = _pos[i];
        int from = position.from;
        bool undirected = position.edge < 0;
        int idx = undirected ? ~position.edge : position.edge;
        const auto& e = _g[from][idx];
        const auto& re = _g[e.to][e.rev];
        if (undirected) {
            return Edge{
                from,
                e.to,
                (e.cap + re.cap) / Cap(2),
                (re.cap - e.cap) / Cap(2)
            };
        }
        return Edge{from, e.to, e.cap + re.cap, re.cap};
    }

    std::vector<Edge> edges() const {
        std::vector<Edge> result;
        result.reserve(_pos.size());
        for (int i = 0; i < int(_pos.size()); i++) result.push_back(get_edge(i));
        return result;
    }

    void change_edge(int i, Cap new_cap, Cap new_flow) {
        assert(0 <= i && i < int(_pos.size()));
        assert(Cap(0) <= new_cap);
        auto& position = _pos[i];
        int from = position.from;
        bool undirected = position.edge < 0;
        int idx = undirected ? ~position.edge : position.edge;
        auto& e = _g[from][idx];
        auto& re = _g[e.to][e.rev];
        if (undirected) {
            assert(new_cap <= std::numeric_limits<Cap>::max() / Cap(2));
            assert(-new_cap <= new_flow && new_flow <= new_cap);
            e.cap = new_cap - new_flow;
            re.cap = new_cap + new_flow;
        } else {
            assert(Cap(0) <= new_flow && new_flow <= new_cap);
            e.cap = new_cap - new_flow;
            re.cap = new_flow;
        }
    }

    Cap max_flow(int s, int t) {
        assert(0 <= s && s < _n);
        assert(0 <= t && t < _n);
        assert(s != t);
        return highest_label_preflow_push(s, t);
    }

    Cap max_flow_push_relabel(int s, int t) {
        assert(0 <= s && s < _n);
        assert(0 <= t && t < _n);
        assert(s != t);
        return highest_label_preflow_push(s, t);
    }

    Cap max_flow_dinic(int s, int t) {
        return max_flow(s, t, std::numeric_limits<Cap>::max());
    }

    Cap max_flow(int s, int t, Cap flow_limit) {
        assert(0 <= s && s < _n);
        assert(0 <= t && t < _n);
        assert(s != t);

        std::vector<int> work(3 * std::size_t(_n));
        int* level = work.data();
        int* iter = level + _n;
        int* queue = iter + _n;
        auto bfs = [&]() -> bool {
            std::fill(level, level + _n, -1);
            int head = 0;
            int tail = 0;
            level[s] = 0;
            queue[tail++] = s;
            while (head != tail) {
                int v = queue[head++];
                for (const auto& e : _g[v]) {
                    if (level[e.to] != -1 || e.cap == Cap(0)) continue;
                    level[e.to] = level[v] + 1;
                    if (e.to == t) return true;
                    queue[tail++] = e.to;
                }
            }
            return level[t] != -1;
        };

        auto dfs = [&](auto&& self, int v, Cap up) -> Cap {
            if (v == s) return up;
            Cap result = Cap(0);
            const int current_level = level[v];
            auto& edges = _g[v];
            const int edge_count = int(edges.size());
            for (int& i = iter[v]; i < edge_count; i++) {
                auto& e = edges[i];
                if (level[e.to] + 1 != current_level) continue;
                auto& reverse = _g[e.to][e.rev];
                if (reverse.cap == Cap(0)) continue;
                Cap d = self(
                    self,
                    e.to,
                    std::min(up - result, reverse.cap)
                );
                if (d == Cap(0)) continue;
                e.cap += d;
                reverse.cap -= d;
                result += d;
                if (result == up) return result;
            }
            level[v] = _n;
            return result;
        };

        Cap flow = 0;
        while (flow < flow_limit && bfs()) {
            std::fill(iter, iter + _n, 0);
            flow += dfs(dfs, t, flow_limit - flow);
        }
        return flow;
    }

    std::vector<bool> min_cut(int s) const {
        assert(0 <= s && s < _n);
        std::vector<bool> visited(_n, false);
        std::vector<int> queue(_n);
        int head = 0;
        int tail = 0;
        visited[s] = true;
        queue[tail++] = s;
        while (head != tail) {
            int v = queue[head++];
            for (const auto& e : _g[v]) {
                if (e.cap == Cap(0) || visited[e.to]) continue;
                visited[e.to] = true;
                queue[tail++] = e.to;
            }
        }
        return visited;
    }
};

}  // namespace flow
}  // namespace m1une


#line 11 "optimization/project_selection.hpp"

namespace m1une {
namespace opt {

template <class T>
struct ProjectSelectionResult {
    bool feasible;
    T max_gain;
    std::vector<bool> selected;

    bool is_feasible() const {
        return feasible;
    }
};

template <class T>
class ProjectSelection {
    static_assert(std::is_integral_v<T> && std::is_signed_v<T>);

    struct Arc {
        int from;
        int to;
        T cap;
    };

    static constexpr int source = -1;
    static constexpr int sink = -2;

    int _project_count;
    int _vertex_count;
    T _offset = T();
    T _finite_cap_sum = T();
    std::vector<Arc> _arcs;
    std::vector<std::pair<int, int>> _hard_arcs;

    void assert_project(int project) const {
        (void)project;
        assert(0 <= project && project < _project_count);
    }

    void assert_vertex(int vertex) const {
        (void)vertex;
        assert(0 <= vertex && vertex < _vertex_count);
    }

    void add_offset(T value) {
        if (value > T()) {
            assert(_offset <= std::numeric_limits<T>::max() - value);
        } else if (value < T()) {
            assert(_offset >= std::numeric_limits<T>::lowest() - value);
        }
        _offset += value;
    }

    T nonnegative_difference(T large, T small) const {
        assert(small <= large);
        if (small < T()) {
            assert(large <= std::numeric_limits<T>::max() + small);
        }
        return large - small;
    }

    void add_arc(int from, int to, T cap) {
        assert(cap >= T());
        if (from == to) return;
        assert(cap <= std::numeric_limits<T>::max() - _finite_cap_sum);
        _finite_cap_sum += cap;
        _arcs.push_back(Arc{from, to, cap});
    }

    void add_hard_arc(int from, int to) {
        if (from == to) return;
        _hard_arcs.emplace_back(from, to);
    }

    void add_vertex_gain(int vertex, T gain_if_selected, T gain_if_unselected) {
        assert_vertex(vertex);
        if (gain_if_selected >= gain_if_unselected) {
            add_offset(gain_if_selected);
            add_arc(source, vertex,
                    nonnegative_difference(gain_if_selected, gain_if_unselected));
        } else {
            add_offset(gain_if_unselected);
            add_arc(vertex, sink,
                    nonnegative_difference(gain_if_unselected, gain_if_selected));
        }
    }

    int add_auxiliary_vertex() {
        return _vertex_count++;
    }

   public:
    ProjectSelection() : ProjectSelection(0) {}

    explicit ProjectSelection(int project_count)
        : _project_count(project_count), _vertex_count(project_count) {
        assert(project_count >= 0);
    }

    int size() const {
        return _project_count;
    }

    void add_gain(int project, T gain_if_selected) {
        add_gain(project, gain_if_selected, T());
    }

    void add_gain(int project, T gain_if_selected, T gain_if_unselected) {
        assert_project(project);
        add_vertex_gain(project, gain_if_selected, gain_if_unselected);
    }

    void add_penalty(int selected_project, int unselected_project, T penalty) {
        assert_project(selected_project);
        assert_project(unselected_project);
        add_arc(selected_project, unselected_project, penalty);
    }

    void add_penalty_if_different(int project_a, int project_b, T penalty) {
        assert_project(project_a);
        assert_project(project_b);
        add_arc(project_a, project_b, penalty);
        add_arc(project_b, project_a, penalty);
    }

    void add_gain_if_same(int project_a, int project_b, T gain) {
        assert(gain >= T());
        add_offset(gain);
        add_penalty_if_different(project_a, project_b, gain);
    }

    void add_hard_implication(int selected_project, int required_project) {
        assert_project(selected_project);
        assert_project(required_project);
        add_hard_arc(selected_project, required_project);
    }

    void force_selected(int project) {
        assert_project(project);
        add_hard_arc(source, project);
    }

    void force_unselected(int project) {
        assert_project(project);
        add_hard_arc(project, sink);
    }

    void add_gain_if_all_selected(const std::vector<int>& projects, T gain) {
        assert(gain >= T());
        for (int project : projects) assert_project(project);
        if (projects.empty()) {
            add_offset(gain);
            return;
        }
        if (projects.size() == 1) {
            add_vertex_gain(projects[0], gain, T());
            return;
        }
        if (projects.size() == 2) {
            add_vertex_gain(projects[0], gain, T());
            add_arc(projects[0], projects[1], gain);
            return;
        }

        int auxiliary = add_auxiliary_vertex();
        add_vertex_gain(auxiliary, gain, T());
        for (int project : projects) add_hard_arc(auxiliary, project);
    }

    void add_gain_if_all_unselected(const std::vector<int>& projects, T gain) {
        assert(gain >= T());
        for (int project : projects) assert_project(project);
        if (projects.empty()) {
            add_offset(gain);
            return;
        }
        if (projects.size() == 1) {
            add_vertex_gain(projects[0], T(), gain);
            return;
        }
        if (projects.size() == 2) {
            add_vertex_gain(projects[0], T(), gain);
            add_arc(projects[1], projects[0], gain);
            return;
        }

        int auxiliary = add_auxiliary_vertex();
        add_vertex_gain(auxiliary, T(), gain);
        for (int project : projects) add_hard_arc(project, auxiliary);
    }

    ProjectSelectionResult<T> solve() const {
        int s = _vertex_count;
        int t = s + 1;
        flow::MaxFlow<T> max_flow(_vertex_count + 2);

        auto vertex_id = [&](int vertex) {
            if (vertex == source) return s;
            if (vertex == sink) return t;
            return vertex;
        };

        for (const auto& arc : _arcs) {
            max_flow.add_edge(vertex_id(arc.from), vertex_id(arc.to), arc.cap);
        }

        T hard_cap = T();
        if (!_hard_arcs.empty()) {
            assert(_finite_cap_sum < std::numeric_limits<T>::max());
            hard_cap = _finite_cap_sum + T(1);
            for (auto [from, to] : _hard_arcs) {
                max_flow.add_edge(vertex_id(from), vertex_id(to), hard_cap);
            }
        }

        T cut_cost =
            _hard_arcs.empty() ? max_flow.max_flow(s, t) : max_flow.max_flow(s, t, hard_cap);
        ProjectSelectionResult<T> result;
        result.feasible = _hard_arcs.empty() || cut_cost < hard_cap;
        result.max_gain = T();
        result.selected.assign(_project_count, false);
        if (!result.feasible) return result;

        assert(_offset >= std::numeric_limits<T>::lowest() + cut_cost);
        result.max_gain = _offset - cut_cost;
        auto source_side = max_flow.min_cut(s);
        for (int project = 0; project < _project_count; project++) {
            result.selected[project] = source_side[project];
        }
        return result;
    }
};

}  // namespace opt
}  // namespace m1une
Back to top page