結果
| 問題 | No.3674 Zero Sum Game |
| コンテスト | |
| ユーザー |
👑  hitonanode
|
| 提出日時 | 2026-10-16 22:51:59 |
| 言語 | C++23 (gcc 15.3.0 + boost 1.92.0 + ACL) |
| 結果 |
AC
不安定
|
| 実行時間 | 105 ms / 3,000 ms |
| + 642µs | |
| コード長 | 95,708 bytes |
| 記録 | |
| コンパイル時間 | 8,539 ms |
| コンパイル使用メモリ | 430,188 KB |
| 実行使用メモリ | 9,928 KB |
| 最終ジャッジ日時 | 2026-10-16 22:52:40 |
| 合計ジャッジ時間 | 17,182 ms |
|
ジャッジサーバーID （参考情報） |
judge2_0 / judge4_0 |
（要ログイン）
| ファイルパターン | 結果 |
|---|---|
| sample | AC * 1 |
| other | AC * 39 |
ソースコード
#include <algorithm>
#include <array>
#include <bitset>
#include <cassert>
#include <chrono>
#include <cmath>
#include <complex>
#include <deque>
#include <forward_list>
#include <fstream>
#include <functional>
#include <iomanip>
#include <ios>
#include <iostream>
#include <limits>
#include <list>
#include <map>
#include <memory>
#include <numeric>
#include <optional>
#include <queue>
#include <random>
#include <set>
#include <sstream>
#include <stack>
#include <string>
#include <tuple>
#include <type_traits>
#include <unordered_map>
#include <unordered_set>
#include <utility>
#include <vector>
using namespace std;
using lint = long long;
using pint = pair<int, int>;
using plint = pair<lint, lint>;
struct fast_ios { fast_ios(){ cin.tie(nullptr), ios::sync_with_stdio(false), cout << fixed << setprecision(20); }; } fast_ios_;
#define ALL(x) (x).begin(), (x).end()
#define FOR(i, begin, end) for(int i=(begin),i##_end_=(end);i<i##_end_;i++)
#define IFOR(i, begin, end) for(int i=(end)-1,i##_begin_=(begin);i>=i##_begin_;i--)
#define REP(i, n) FOR(i,0,n)
#define IREP(i, n) IFOR(i,0,n)
template <typename T> bool chmax(T &m, const T q) { return m < q ? (m = q, true) : false; }
template <typename T> bool chmin(T &m, const T q) { return m > q ? (m = q, true) : false; }
const std::vector<std::pair<int, int>> grid_dxs{{1, 0}, {-1, 0}, {0, 1}, {0, -1}};
int floor_lg(long long x) { return x <= 0 ? -1 : 63 - __builtin_clzll(x); }
template <class T1, class T2> T1 floor_div(T1 num, T2 den) { return (num > 0 ? num / den : -((-num + den - 1) / den)); }
template <class T1, class T2> std::pair<T1, T2> operator+(const std::pair<T1, T2> &l, const std::pair<T1, T2> &r) { return std::make_pair(l.first + r.first, l.second + r.second); }
template <class T1, class T2> std::pair<T1, T2> operator-(const std::pair<T1, T2> &l, const std::pair<T1, T2> &r) { return std::make_pair(l.first - r.first, l.second - r.second); }
template <class T> std::vector<T> sort_unique(std::vector<T> vec) { sort(vec.begin(), vec.end()), vec.erase(unique(vec.begin(), vec.end()), vec.end()); return vec; }
template <class T> int arglb(const std::vector<T> &v, const T &x) { return std::distance(v.begin(), std::lower_bound(v.begin(), v.end(), x)); }
template <class T> int argub(const std::vector<T> &v, const T &x) { return std::distance(v.begin(), std::upper_bound(v.begin(), v.end(), x)); }
template <class IStream, class T> IStream &operator>>(IStream &is, std::vector<T> &vec) { for (auto &v : vec) is >> v; return is; }
template <class OStream, class T> OStream &operator<<(OStream &os, const std::vector<T> &vec);
template <class OStream, class T, size_t sz> OStream &operator<<(OStream &os, const std::array<T, sz> &arr);
template <class OStream, class T, class TH> OStream &operator<<(OStream &os, const std::unordered_set<T, TH> &vec);
template <class OStream, class T, class U> OStream &operator<<(OStream &os, const pair<T, U> &pa);
template <class OStream, class T> OStream &operator<<(OStream &os, const std::deque<T> &vec);
template <class OStream, class T> OStream &operator<<(OStream &os, const std::set<T> &vec);
template <class OStream, class T> OStream &operator<<(OStream &os, const std::multiset<T> &vec);
template <class OStream, class T> OStream &operator<<(OStream &os, const std::unordered_multiset<T> &vec);
template <class OStream, class T, class U> OStream &operator<<(OStream &os, const std::pair<T, U> &pa);
template <class OStream, class TK, class TV> OStream &operator<<(OStream &os, const std::map<TK, TV> &mp);
template <class OStream, class TK, class TV, class TH> OStream &operator<<(OStream &os, const std::unordered_map<TK, TV, TH> &mp);
template <class OStream, class... T> OStream &operator<<(OStream &os, const std::tuple<T...> &tpl);
template <class OStream, class T> OStream &operator<<(OStream &os, const std::vector<T> &vec) { os << '['; for (auto v : vec) os << v << ','; os << ']'; return os; }
template <class OStream, class T, size_t sz> OStream &operator<<(OStream &os, const std::array<T, sz> &arr) { os << '['; for (auto v : arr) os << v << ','; os << ']'; return os; }
template <class... T> std::istream &operator>>(std::istream &is, std::tuple<T...> &tpl) { std::apply([&is](auto &&... args) { ((is >> args), ...);}, tpl); return is; }
template <class OStream, class... T> OStream &operator<<(OStream &os, const std::tuple<T...> &tpl) { os << '('; std::apply([&os](auto &&... args) { ((os << args << ','), ...);}, tpl); return os << ')'; }
template <class OStream, class T, class TH> OStream &operator<<(OStream &os, const std::unordered_set<T, TH> &vec) { os << '{'; for (auto v : vec) os << v << ','; os << '}'; return os; }
template <class OStream, class T> OStream &operator<<(OStream &os, const std::deque<T> &vec) { os << "deq["; for (auto v : vec) os << v << ','; os << ']'; return os; }
template <class OStream, class T> OStream &operator<<(OStream &os, const std::set<T> &vec) { os << '{'; for (auto v : vec) os << v << ','; os << '}'; return os; }
template <class OStream, class T> OStream &operator<<(OStream &os, const std::multiset<T> &vec) { os << '{'; for (auto v : vec) os << v << ','; os << '}'; return os; }
template <class OStream, class T> OStream &operator<<(OStream &os, const std::unordered_multiset<T> &vec) { os << '{'; for (auto v : vec) os << v << ','; os << '}'; return os; }
template <class OStream, class T, class U> OStream &operator<<(OStream &os, const std::pair<T, U> &pa) { return os << '(' << pa.first << ',' << pa.second << ')'; }
template <class OStream, class TK, class TV> OStream &operator<<(OStream &os, const std::map<TK, TV> &mp) { os << '{'; for (auto v : mp) os << v.first << "=>" << v.second << ','; os << '}'; return os; }
template <class OStream, class TK, class TV, class TH> OStream &operator<<(OStream &os, const std::unordered_map<TK, TV, TH> &mp) { os << '{'; for (auto v : mp) os << v.first << "=>" << v.second << ','; os << '}'; return os; }
#ifdef HITONANODE_LOCAL
const string COLOR_RESET = "\033[0m", BRIGHT_CYAN = "\033[1;36m", NORMAL_FAINT = "\033[0;2m";
#define dbg(x) std::cerr << BRIGHT_CYAN << #x << COLOR_RESET << " = " << (x) << NORMAL_FAINT << " (L" << __LINE__ << ") " << __FILE__ << COLOR_RESET << std::endl
#define dbgif(cond, x) ((cond) ? std::cerr << BRIGHT_CYAN << #x << COLOR_RESET << " = " << (x) << NORMAL_FAINT << " (L" << __LINE__ << ") " << __FILE__ << COLOR_RESET << std::endl : std::cerr)
#else
#define dbg(x) ((void)0)
#define dbgif(cond, x) ((void)0)
#endif
#ifndef mip_HPP
#define mip_HPP
#include <algorithm>
#include <chrono>
#include <cmath>
#include <cstdint>
#include <functional>
#include <limits>
#include <memory>
#include <numeric>
#include <stdexcept>
#include <utility>
#include <vector>
namespace cpopt {
inline constexpr double inf = std::numeric_limits<double>::infinity();
enum class Status { Optimal, Infeasible, Unbounded, Limit, NumericalError, RelaxationUnbounded };
enum class Sense { Minimize, Maximize };
enum class VarType { Continuous, Integer, Binary };
enum class Relation { LessEqual, Equal, GreaterEqual };
using Terms = std::vector<std::pair<int, double>>;
inline const char *status_name(Status s) {
switch (s) {
case Status::Optimal: return "optimal";
case Status::Infeasible: return "infeasible";
case Status::Unbounded: return "unbounded";
case Status::Limit: return "limit";
case Status::NumericalError: return "numerical_error";
case Status::RelaxationUnbounded: return "relaxation_unbounded";
}
return "unknown";
}
struct Options {
double time_limit = inf; // 秒。前処理を含む solve() ごとの経過時間。
std::uint64_t node_limit = std::numeric_limits<std::uint64_t>::max();
std::uint64_t iteration_limit = std::numeric_limits<std::uint64_t>::max();
double feasibility_tolerance = 1e-7;
double integrality_tolerance = 1e-7;
double pivot_tolerance = 1e-11;
double absolute_gap = 1e-7;
double relative_gap = 0;
bool presolve = true;
bool scaling = true;
bool dualize = true; // 連続 LP で有利な形なら双対問題を解く。
bool warm_start = true;
bool reuse_root = true; // 同じ Model の再求解で、最適な根基底を保存・再利用。
bool pseudocost = true;
bool rounding_heuristic = true;
bool integer_objective = true;
bool reduced_cost_fixing = true; // 暫定解を改善できない非基底整数変数を固定・列削除。
int root_cut_rounds = -1; // -1: 小数の整数変数が1個のときだけ3巡、0: 無効。
int cuts_per_round = 8;
int strong_branching_candidates = 4;
int strong_branching_depth = 3;
int strong_branching_iterations = 40;
std::vector<double> initial_solution; // 元の変数順。実行可能性を検査してから使う。
};
struct Statistics {
std::uint64_t nodes = 0, pivots = 0, strong_branching_pivots = 0;
std::uint64_t heuristic_solutions = 0;
int cuts = 0;
int reduced_cost_fixings = 0;
int numerical_restarts = 0;
bool solved_dual = false;
bool reused_root = false;
std::uint64_t root_pivots = 0; // 根の初回求解・再最適化・復旧の総 pivot 数。
int original_variables = 0, original_constraints = 0;
int presolved_variables = 0, presolved_constraints = 0;
double seconds = 0;
};
struct Result {
Status status = Status::NumericalError;
bool has_solution = false;
double objective = std::numeric_limits<double>::quiet_NaN();
// 最大化なら最適値の上界、最小化なら下界。解がなくても有効。
double best_bound = std::numeric_limits<double>::quiet_NaN();
std::vector<double> x;
Statistics stats;
double gap() const { return has_solution ? std::abs(best_bound - objective) : inf; }
};
namespace detail {
// 部分的に構築した前処理結果を誤って利用せず、保持済みの暫定解へ戻る。
struct TimeLimitReached {};
struct Variable {
double lower, upper;
VarType type;
};
struct Row {
Terms terms;
double rhs;
}; // 常に <=
struct Budget {
const Options &opt;
Statistics &stats;
std::chrono::steady_clock::time_point start = std::chrono::steady_clock::now();
double elapsed() const {
return std::chrono::duration<double>(std::chrono::steady_clock::now() - start).count();
}
static constexpr std::size_t check_interval = 4096;
const bool timed = std::isfinite(opt.time_limit);
std::size_t work_left = 0;
// 係数ごとには時計を読まない。処理量を積算し、約4096要素ごとに読む。
void checkpoint(std::size_t work = check_interval) {
if (!timed) return;
if (work < work_left) {
work_left -= work;
return;
}
work_left = check_interval;
if (elapsed() >= opt.time_limit) throw TimeLimitReached{};
}
// 密なループはブロック内に条件分岐を入れず、自動ベクトル化を維持する。
template <class F> void chunks(std::size_t size, F f) {
if (!timed) {
f(0, size);
return;
}
for (std::size_t begin = 0; begin < size; begin += check_interval) {
const auto end = std::min(size, begin + check_interval);
checkpoint(end - begin);
f(begin, end);
}
}
template <class F> static bool scan(Budget *budget, std::size_t size, F f) {
if (!budget || !budget->timed) return f(0, size);
for (std::size_t begin = 0; begin < size; begin += check_interval) {
const auto end = std::min(size, begin + check_interval);
budget->checkpoint(end - begin);
if (!f(begin, end)) return false;
}
return true;
}
template <class T> std::vector<T> copy(const std::vector<T> &source) {
if (!timed) return source;
checkpoint(1);
std::vector<T> result;
result.reserve(source.size());
chunks(source.size(), [&](std::size_t begin, std::size_t end) {
result.insert(result.end(), source.begin() + begin, source.begin() + end);
});
return result;
}
template <class T> std::vector<T> filled(std::size_t size, const T &value = T{}) {
if (!timed) return std::vector<T>(size, value);
checkpoint(1);
std::vector<T> result;
result.reserve(size);
chunks(size, [&](std::size_t begin, std::size_t end) {
result.insert(result.end(), end - begin, value);
});
return result;
}
bool stop() const {
return stats.pivots >= opt.iteration_limit ||
(std::isfinite(opt.time_limit) && elapsed() >= opt.time_limit);
}
};
// 辞書: x_B = rhs - D x_N, z = rhs_z - D_z x_N。
// 初期列数は元変数数 + 人工変数 + RHS。固定列は後から削除できる。
// slack の列を別に持たず交換する。
struct Tableau {
int n = 0, m = 0, stride = 0, next_id = 0, original_n = 0;
static constexpr int fixed_zero = std::numeric_limits<int>::max();
bool primal_feasible = false;
double objective_scale = 1;
std::vector<double> dense;
std::vector<std::shared_ptr<std::vector<double>>> data;
std::vector<int> basic, nonbasic;
std::vector<int> location; // fixed_zero: 0固定、それ以外は非負: 列、負: -行-1。人工変数を除く。
std::vector<int> pivot_columns; // 疎な pivot 行の非ゼロ列。作業領域なので複製不要。
std::vector<double *> pivot_rows; // detach 済みの行。pivot 中だけ有効。
const double *row(int i) const {
return data.empty() ? dense.data() + std::size_t(i) * stride : data[i]->data();
}
double *write_row(int i, Budget *budget = nullptr) {
if (data.empty()) return dense.data() + std::size_t(i) * stride;
if (data[i].use_count() != 1)
data[i] =
std::make_shared<std::vector<double>>(budget ? budget->copy(*data[i]) : *data[i]);
return data[i]->data();
}
double value() const { return row(m)[n + 1] / objective_scale; }
void scale_objective(const std::vector<double> &c) {
double largest = 0;
for (double v : c) largest = std::max(largest, std::abs(v));
int e = 0;
if (largest) std::frexp(largest, &e);
// 小さすぎる目的係数だけ拡大する。大きい目的を縮めると、例えば
// 1e12*x + y の y の寄与が pivot 許容誤差より小さくなってしまう。
objective_scale = std::ldexp(1.0, -std::max(-1000, std::min(0, e)));
}
Tableau() = default;
Tableau(const Tableau &other) { *this = other; }
Tableau &operator=(const Tableau &other) {
if (this == &other) return *this;
n = other.n;
m = other.m;
stride = other.stride;
next_id = other.next_id;
original_n = other.original_n;
objective_scale = other.objective_scale;
primal_feasible = other.primal_feasible;
// 複製直後の分枝行1本で再確保しない。値のコピーと行追加で二重に
// 配列をコピーするコストを避ける（探索・分枝の順序は同じ）。
if (!other.data.empty()) data.reserve(other.data.size() + 1);
data = other.data;
dense.reserve(other.dense.size() + (other.dense.empty() ? 0 : std::size_t(stride) * 2));
dense.assign(other.dense.begin(), other.dense.end());
basic.reserve(other.basic.size() + 1);
location.reserve(other.location.size() + 1);
basic = other.basic;
nonbasic = other.nonbasic;
location = other.location;
return *this;
}
// 再求解の保存・復元時は、既存の不変スナップショットを壊さず複製する。
Tableau clone(Budget &budget) const {
Tableau t;
t.n = n;
t.m = m;
t.stride = stride;
t.next_id = next_id;
t.original_n = original_n;
t.primal_feasible = primal_feasible;
t.objective_scale = objective_scale;
budget.checkpoint();
t.data = budget.copy(data);
t.dense.reserve(dense.size() + (dense.empty() ? 0 : std::size_t(stride) * 2));
budget.chunks(dense.size(), [&](std::size_t begin, std::size_t end) {
t.dense.insert(t.dense.end(), dense.begin() + begin, dense.begin() + end);
});
t.basic = budget.copy(basic);
t.nonbasic = budget.copy(nonbasic);
t.location = budget.copy(location);
return t;
}
Tableau(Tableau &&) = default;
Tableau &operator=(Tableau &&) = default;
static void copy_coefficients(const std::vector<double> &a, double *dest, Budget *budget) {
auto copy = [&](std::size_t begin, std::size_t end) {
std::copy(a.begin() + begin, a.begin() + end, dest + begin);
};
if (budget)
budget->chunks(a.size(), copy);
else
copy(0, a.size());
}
static void copy_coefficients(const Terms &a, double *dest, Budget *budget) {
auto copy = [&](std::size_t begin, std::size_t end) {
for (auto k = begin; k < end; ++k) dest[a[k].first] = a[k].second;
};
if (budget)
budget->chunks(a.size(), copy);
else
copy(0, a.size());
}
template <class Matrix>
Tableau(const Matrix &a, const std::vector<double> &b, const std::vector<double> &c,
Budget *budget = nullptr)
: n(int(c.size())), m(int(b.size())), stride(n + 2), next_id(n + m), original_n(n),
basic(m), nonbasic(n + 1), location(n + m) {
if (stride < 128)
dense = budget ? budget->filled<double>(std::size_t(m + 2) * stride)
: std::vector<double>(std::size_t(m + 2) * stride);
else {
data.reserve(m + 2);
for (int i = 0; i < m + 2; ++i)
data.push_back(std::make_shared<std::vector<double>>(
budget ? budget->filled<double>(stride) : std::vector<double>(stride)));
}
for (int i = 0; i < m; ++i) {
double *q = write_row(i);
copy_coefficients(a[i], q, budget);
q[n] = -1;
q[n + 1] = b[i];
basic[i] = n + i;
location[n + i] = -i - 1;
}
scale_objective(c);
double *z = write_row(m);
for (int j = 0; j < n; ++j) {
nonbasic[j] = j;
location[j] = j;
z[j] = -c[j] * objective_scale;
}
nonbasic[n] = -1;
write_row(m + 1)[n] = 1;
}
void pivot(int r, int s, Budget &budget) {
double *p = write_row(r, &budget);
pivot_rows.clear();
if (!data.empty())
for (int i = 0; i < m + 2; ++i)
if (i != r && row(i)[s] != 0) { pivot_rows.push_back(write_row(i, &budget)); }
const double inv = 1 / p[s];
// フロー等では辞書を密に保存していても pivot 行は疎なことが多い。
// 密な行は早めに走査を打ち切り、従来の連続ループを使う。
pivot_columns.clear();
const int sparse_limit = stride >= 64 ? stride / 4 : 0;
if (sparse_limit) {
for (int j = 0; j < stride; ++j)
if (p[j] != 0) {
pivot_columns.push_back(j);
if (int(pivot_columns.size()) > sparse_limit) break;
}
}
const bool sparse = sparse_limit && int(pivot_columns.size()) <= sparse_limit;
// 連続メモリでの rank-one 更新。ゼロ行は丸ごと省く。
auto update = [&](double *q) {
const double f = q[s] * inv;
if (f == 0) return;
if (sparse) {
for (int j : pivot_columns) q[j] -= f * p[j];
} else {
for (int j = 0; j < stride; ++j) q[j] -= f * p[j];
}
q[s] = -f;
};
// 格納方式の分岐と共有ポインタの参照は、数値更新の内側に持ち込まない。
if (data.empty()) {
for (double *q = dense.data(); q != p; q += stride) update(q);
double *end = dense.data() + dense.size();
for (double *q = p + stride; q != end; q += stride) update(q);
} else
for (double *q : pivot_rows) update(q);
for (int j = 0; j < stride; ++j) p[j] *= inv;
p[s] = inv;
std::swap(basic[r], nonbasic[s]);
if (basic[r] >= 0) location[basic[r]] = -r - 1;
if (nonbasic[s] >= 0) location[nonbasic[s]] = s;
++budget.stats.pivots;
}
Status primal(int objective_row, Budget &budget) {
const double eps = budget.opt.pivot_tolerance;
int degenerate = 0;
bool bland = false;
while (true) {
if (budget.stop()) return Status::Limit;
const double *z = row(objective_row);
int s = -1;
for (int j = 0; j <= n; ++j) {
if (nonbasic[j] == -1 || z[j] >= -eps) continue;
if (s < 0 || (bland ? nonbasic[j] < nonbasic[s] : z[j] < z[s])) s = j;
}
if (s < 0) return Status::Optimal;
int r = -1;
double ratio = inf;
for (int i = 0; i < m; ++i) {
const double *q = row(i);
if (q[s] <= eps) continue;
const double v = std::max(0.0, q[n + 1]) / q[s];
if (v < ratio - eps ||
(std::abs(v - ratio) <= eps && (r < 0 || basic[i] < basic[r]))) {
r = i;
ratio = v;
}
}
if (r < 0) return Status::Unbounded;
if (ratio <= eps)
++degenerate;
else {
degenerate = 0;
bland = false;
}
if (degenerate > 32) bland = true; // 退化が続いたら循環を避ける。
pivot(r, s, budget);
if (!std::isfinite(value())) return Status::NumericalError;
}
}
Status solve(Budget &budget) {
if (budget.stop()) return Status::Limit;
int r = -1;
for (int i = 0; i < m; ++i)
if (r < 0 || row(i)[n + 1] < row(r)[n + 1]) r = i;
if (r >= 0 && row(r)[n + 1] < -budget.opt.pivot_tolerance) {
// 最小費用フローや被覆では、slack 基底が既に双対実行可能。
// 人工変数で全行を混ぜる Phase I を省き、直接実行可能性を回復する。
bool dual_feasible = true;
for (int j = 0; j <= n; ++j)
if (nonbasic[j] != -1 && row(m)[j] < 0) dual_feasible = false;
if (dual_feasible) return dual(budget);
pivot(r, n, budget);
Status s = primal(m + 1, budget);
if (s == Status::Limit || s == Status::NumericalError) return s;
if (s != Status::Optimal) return Status::NumericalError;
if (row(m + 1)[n + 1] < -budget.opt.pivot_tolerance) return Status::Infeasible;
for (int i = 0; i < m; ++i)
if (basic[i] == -1) {
int jbest = -1;
for (int j = 0; j <= n; ++j)
if (nonbasic[j] != -1 && std::abs(row(i)[j]) > budget.opt.pivot_tolerance &&
(jbest < 0 || std::abs(row(i)[j]) > std::abs(row(i)[jbest])))
jbest = j;
if (jbest >= 0) pivot(i, jbest, budget);
}
}
// 人工変数が非基底なら選ばない。冗長行の基底に残った場合も値は 0。
primal_feasible = true;
return primal(m, budget);
}
std::vector<double> solution() const {
std::vector<double> x(original_n);
for (int j = 0; j < original_n; ++j)
if (location[j] < 0) x[j] = row(-location[j] - 1)[n + 1];
return x;
}
// 正の reduced cost と整数の最小移動幅で、改善解では0と証明できる列を除く。
// 元変数/slackのIDは保ち、辞書の物理的な列数だけを縮める。
int fix_nonbasic(const std::vector<double> &step, double gap, Budget &budget,
int minimum_removed = 1) {
// 多くのノードでは削除しない。pivot の作業領域に候補だけを記録し、
// 閾値に達してから location と残す列の配列を作る。
pivot_columns.clear();
const double margin = budget.opt.feasibility_tolerance * (1 + std::abs(value()));
const double *z = row(m);
for (int j = 0; j <= n; ++j) {
budget.checkpoint(1);
const int id = nonbasic[j];
if (id >= 0 && step[id] > 0 && z[j] / objective_scale * step[id] > gap + margin)
pivot_columns.push_back(j);
}
const int removed = int(pivot_columns.size());
if (removed < minimum_removed) return 0;
auto positions = budget.copy(location);
std::vector<int> keep;
keep.reserve(n + 1 - removed);
int k = 0;
for (int j = 0; j <= n; ++j) {
budget.checkpoint(1);
if (k < removed && pivot_columns[k] == j) {
positions[nonbasic[j]] = fixed_zero;
++k;
} else
keep.push_back(j);
}
const int width = int(keep.size()) + 1;
std::vector<std::shared_ptr<std::vector<double>>> compact;
std::vector<double> compact_dense;
if (width < 128)
compact_dense = budget.filled<double>(std::size_t(m + 2) * width);
else
compact.reserve(m + 2);
for (int i = 0; i < m + 2; ++i) {
auto values = width < 128 ? std::vector<double>{} : budget.filled<double>(width);
double *dest =
width < 128 ? compact_dense.data() + std::size_t(i) * width : values.data();
const double *source = row(i);
budget.chunks(keep.size(), [&](std::size_t begin, std::size_t end) {
for (auto j = begin; j < end; ++j) dest[j] = source[keep[j]];
});
dest[width - 1] = source[n + 1];
if (width >= 128)
compact.push_back(std::make_shared<std::vector<double>>(std::move(values)));
}
std::vector<int> columns;
columns.reserve(keep.size());
for (int j : keep) {
const int id = nonbasic[j];
if (id >= 0) positions[id] = int(columns.size());
columns.push_back(id);
}
data = std::move(compact);
dense = std::move(compact_dense);
nonbasic = std::move(columns);
location = std::move(positions);
n -= removed;
stride = width;
return removed;
}
// 最適辞書に a*x <= rhs を追加。目的行の双対実行可能性は保存される。
template <class Coefficients = Terms>
void add_constraint(const Coefficients &a, double rhs, Budget *budget = nullptr,
double multiplier = 1) {
primal_feasible = false;
// 時間割込みのない狭い辞書では、追加先に直接組み立てて一時行を省く。
// budget ありの場合は従来通り、完成した行だけを辞書へ追加する。
const bool inplace = data.empty() && !budget;
std::vector<double> buffer;
double *q;
if (inplace) {
dense.resize(std::size_t(m + 3) * stride);
std::copy_backward(write_row(m), write_row(m + 2), write_row(m + 3));
q = write_row(m);
std::fill(q, q + stride, 0);
} else {
buffer.resize(stride);
q = buffer.data();
}
q[n + 1] = rhs;
for (const auto &term : a) {
if (budget) budget->checkpoint(1);
const int var = term.first;
const double coef = term.second * multiplier;
const int pos = location[var];
if (pos == fixed_zero) continue;
if (pos >= 0)
q[pos] += coef;
else {
const double *p = row(-pos - 1);
auto subtract = [&](std::size_t begin, std::size_t end) {
for (auto j = begin; j < end; ++j) q[j] -= coef * p[j];
};
if (budget)
budget->chunks(stride, subtract);
else
subtract(0, stride);
}
}
if (data.empty() && !inplace) {
dense.resize(std::size_t(m + 3) * stride);
std::copy_backward(write_row(m), write_row(m + 2), write_row(m + 3));
std::copy(buffer.begin(), buffer.end(), write_row(m));
} else if (!data.empty())
data.insert(data.begin() + m, std::make_shared<std::vector<double>>(std::move(buffer)));
basic.push_back(next_id++);
location.push_back(-m - 1);
++m;
}
Status dual(Budget &budget, int max_pivots = std::numeric_limits<int>::max()) {
const double eps = budget.opt.pivot_tolerance;
int iter = 0;
int degenerate = 0;
bool bland = false;
while (true) {
if (budget.stop()) return Status::Limit;
int r = -1;
for (int i = 0; i < m; ++i)
if (row(i)[n + 1] < -eps &&
(r < 0 || (bland ? basic[i] < basic[r] : row(i)[n + 1] < row(r)[n + 1])))
r = i;
if (r < 0) {
primal_feasible = true;
return Status::Optimal;
}
if (iter++ >= max_pivots) return Status::Limit;
int s = -1;
double ratio = inf;
const double *q = row(r);
const double *z = row(m);
for (int j = 0; j <= n; ++j) {
if (nonbasic[j] == -1 || q[j] >= -eps) continue;
const double v = std::max(0.0, z[j]) / -q[j];
if (v < ratio - eps ||
(std::abs(v - ratio) <= eps && (s < 0 || nonbasic[j] < nonbasic[s]))) {
s = j;
ratio = v;
}
}
if (s < 0) return Status::Infeasible;
if (ratio <= eps)
++degenerate;
else {
degenerate = 0;
bland = false;
}
if (degenerate > 32) bland = true;
pivot(r, s, budget);
if (!std::isfinite(value())) return Status::NumericalError;
}
}
};
// 元変数は固定・平行移動・反転・自由変数の分割のいずれかで、対応は高々2項。
// 可変長配列を変数ごとに確保せず、根の保存・復元でもまとめてコピーする。
struct Mapping {
struct Columns {
std::pair<int, double> entries[2];
int count = 0;
auto begin() const { return entries; }
auto end() const { return entries + count; }
int size() const { return count; }
bool empty() const { return count == 0; }
const auto &operator[](int j) const { return entries[j]; }
void emplace_back(int j, double a) { entries[count++] = {j, a}; }
} terms;
double offset = 0;
};
// 変換済みの係数行は追加後に変更しない。根の保存・復元では行を共有し、
// 追加行・カットのリストだけを独立に持つ。Model のコピーにも変更が漏れない。
struct ImmutableRows {
std::vector<std::shared_ptr<const Terms>> rows;
std::size_t size() const { return rows.size(); }
void reserve(std::size_t n) { rows.reserve(n); }
const Terms &operator[](std::size_t i) const { return *rows[i]; }
void push_terms(Terms row) { rows.push_back(std::make_shared<const Terms>(std::move(row))); }
void push_back(const std::vector<double> &row) {
Terms terms;
for (int j = 0; j < int(row.size()); ++j)
if (row[j]) terms.emplace_back(j, row[j]);
push_terms(std::move(terms));
}
};
struct Prepared {
std::vector<Variable> vars;
std::vector<Row> rows;
std::vector<Mapping> map;
ImmutableRows a;
std::vector<double> b, c;
double offset = 0;
bool infeasible = false;
Prepared clone(Budget &budget) const {
Prepared p;
p.vars = budget.copy(vars);
p.map = budget.copy(map);
p.b = budget.copy(b);
p.c = budget.copy(c);
p.offset = offset;
p.infeasible = infeasible;
p.a.rows = budget.copy(a.rows);
// rows は前処理だけで使う作業領域。保存・再求解には不要。
return p;
}
};
// Gomory mixed-integer cut。教科書の式から生成し、非基底 slack を元の行で
// 消去して保持する。これにより cold start でも同じカットを再現できる。
inline int gomory_cuts(Tableau &tab, Prepared &p, const std::vector<bool> &integer, Budget &budget) {
const auto &opt = budget.opt;
struct Cut {
Terms terms;
double score;
};
std::vector<Cut> cuts;
for (int i = 0; i < tab.m; ++i) {
budget.checkpoint(1);
int var = tab.basic[i];
if (var < 0 || var >= tab.original_n || !integer[var]) continue;
const double *row = tab.row(i);
const double f = row[tab.n + 1] - std::floor(row[tab.n + 1]);
if (f < .01 || f > .99) continue;
Terms terms;
double norm = 0, largest = 0;
for (int j = 0; j <= tab.n; ++j) {
budget.checkpoint(1);
const int id = tab.nonbasic[j];
if (id < 0) continue;
const double a = row[j];
double alpha;
if (id < tab.original_n && integer[id]) {
const double fj = a - std::floor(a);
alpha = fj <= f ? fj / f : (1 - fj) / (1 - f);
} else
alpha = a >= 0 ? a / f : -a / (1 - f);
if (alpha == 0) continue;
terms.emplace_back(id, -alpha);
norm += alpha * alpha;
largest = std::max(largest, alpha);
}
// 極端な係数の行は切らず、通常の分枝で処理する。
if (!std::isfinite(norm) || largest > 1e6) continue;
cuts.push_back({std::move(terms), 1 / std::max(1e-12, std::sqrt(norm))});
}
std::stable_sort(
cuts.begin(), cuts.end(), [](const Cut &a, const Cut &b) { return a.score > b.score; });
int added = 0;
for (const auto &cut : cuts) {
budget.checkpoint();
if (added >= opt.cuts_per_round) break;
std::vector<long double> a(tab.original_n);
// 丸め誤差に対してわずかに外側へ緩める。
long double rhs = -1 + 1e-8;
for (const auto &term : cut.terms) {
const int id = term.first;
const double coef = term.second;
if (id < tab.original_n)
a[id] += coef;
else {
const int r = id - tab.original_n;
rhs -= (long double)coef * p.b[r];
budget.chunks(p.a[r].size(), [&](std::size_t begin, std::size_t end) {
for (auto k = begin; k < end; ++k) {
auto [j, v] = p.a[r][k];
a[j] -= (long double)coef * v;
}
});
}
}
Terms coefficients;
double max_abs = 0;
budget.chunks(a.size(), [&](std::size_t begin, std::size_t end) {
for (auto j = begin; j < end; ++j) {
const double v = double(a[j]);
max_abs = std::max(max_abs, std::abs(v));
if (v) coefficients.emplace_back(int(j), v);
}
});
if (!std::isfinite(double(rhs)) || !std::isfinite(max_abs) || max_abs > 1e8) continue;
p.a.push_terms(std::move(coefficients));
p.b.push_back(double(rhs));
tab.add_constraint(cut.terms, -1 + 1e-8);
++added;
}
return added;
}
inline double bound_round(double v, bool lower, double tol) {
return lower ? std::ceil(v - tol) : std::floor(v + tol);
}
// 各行の区間活動量から境界を伝播する。無限境界も扱う。
inline bool propagate(const std::vector<Row> &rows, std::vector<Variable> &vars, Budget &budget,
int rounds = 8) {
const auto &opt = budget.opt;
for (int pass = 0; pass < rounds; ++pass) {
bool changed = false;
for (const Row &r : rows) {
budget.checkpoint(1);
long double sum = 0;
int infinite = 0;
budget.chunks(r.terms.size(), [&](std::size_t begin, std::size_t end) {
for (auto k = begin; k < end; ++k) {
const auto [j, a] = r.terms[k];
const double b = a > 0 ? vars[j].lower : vars[j].upper;
if (std::isfinite(b))
sum += (long double)a * b;
else
++infinite;
}
});
const long double slack =
opt.feasibility_tolerance + 1e-12L * (std::abs(sum) + std::abs(r.rhs));
if (!infinite && sum > (long double)r.rhs + slack) return false;
if (!Budget::scan(&budget, r.terms.size(), [&](std::size_t begin, std::size_t end) {
for (auto k = begin; k < end; ++k) {
const auto [j, a] = r.terms[k];
Variable &v = vars[j];
const double old = a > 0 ? v.lower : v.upper;
const bool finite = std::isfinite(old);
if (infinite - !finite != 0) continue;
const long double other = sum - (finite ? (long double)a * old : 0);
double b = double(((long double)r.rhs - other) / a);
if (!std::isfinite(b)) continue;
const bool lower = a < 0;
if (v.type != VarType::Continuous)
b = bound_round(b, lower, opt.integrality_tolerance);
double &target = lower ? v.lower : v.upper;
if (lower ? b > target : b < target) {
// 各行内での活動量は更新前の境界に対応する。古い境界からも有効な推論。
// 微小な連続境界の更新による収束ループを抑える。
if (std::isfinite(target) &&
std::abs(b - target) <= 1e-10 * (1 + std::abs(b)))
continue;
target = b;
changed = true;
}
if (v.lower > v.upper + opt.feasibility_tolerance) return false;
}
return true;
}))
return false;
}
if (!changed) break;
}
return true;
}
} // namespace detail
class Model {
std::vector<detail::Variable> vars_;
std::vector<detail::Row> rows_;
std::vector<double> objective_;
Sense sense_ = Sense::Minimize;
double constant_ = 0;
struct RootCache {
detail::Prepared prepared;
detail::Tableau tableau;
std::size_t rows = 0;
bool presolve = false, scaling = false, dualize = false, relax = false;
double feasibility = 0, integrality = 0, pivot = 0;
bool matches(const Options &o, bool r) const {
return presolve == o.presolve && scaling == o.scaling && dualize == o.dualize &&
relax == r && feasibility == o.feasibility_tolerance &&
integrality == o.integrality_tolerance && pivot == o.pivot_tolerance;
}
};
// Model のコピー間でも読み取り専用。追加行やカットは必ず別の作業用辞書に入れる。
mutable std::shared_ptr<const RootCache> root_cache_;
void check_var(int j) const {
if (j < 0 || j >= int(vars_.size())) throw std::invalid_argument("invalid variable index");
}
Terms normalize(Terms terms) const {
for (auto [j, a] : terms) {
check_var(j);
if (!std::isfinite(a)) throw std::invalid_argument("non-finite coefficient");
}
std::sort(terms.begin(), terms.end());
std::size_t size = 0;
for (auto [j, a] : terms) {
if (size && terms[size - 1].first == j)
terms[size - 1].second += a;
else
terms[size++] = {j, a};
}
terms.resize(size);
terms.erase(std::remove_if(terms.begin(), terms.end(), [](auto p) { return p.second == 0; }),
terms.end());
for (auto p : terms)
if (!std::isfinite(p.second)) throw std::invalid_argument("coefficient overflow");
return terms;
}
detail::Prepared prepare(detail::Budget &budget, bool relax) const {
const auto &opt = budget.opt;
detail::Prepared p;
p.vars = budget.copy(vars_);
if (!budget.timed)
p.rows = rows_;
else {
p.rows.reserve(rows_.size());
for (const auto &r : rows_) {
budget.checkpoint(1);
p.rows.push_back({budget.copy(r.terms), r.rhs});
}
}
for (auto &v : p.vars) {
budget.checkpoint(1);
if (relax) v.type = VarType::Continuous;
if (v.type != VarType::Continuous) {
v.lower = detail::bound_round(v.lower, true, opt.integrality_tolerance);
v.upper = detail::bound_round(v.upper, false, opt.integrality_tolerance);
}
if (v.lower > v.upper) {
p.infeasible = true;
return p;
}
}
if (opt.presolve) {
// 整数係数だけの行は gcd で RHS を切り下げられる。
for (auto &r : p.rows) {
budget.checkpoint(1);
long long g = 0;
bool integral = true;
for (auto [j, a] : r.terms) {
budget.checkpoint(1);
if (p.vars[j].type == VarType::Continuous || std::abs(a) > 1e9 ||
a != std::round(a)) {
integral = false;
break;
}
g = std::gcd(g, (long long)std::llround(std::abs(a)));
}
if (integral && g > 0) {
r.rhs = std::floor(r.rhs / g + opt.integrality_tolerance);
budget.chunks(r.terms.size(), [&](std::size_t begin, std::size_t end) {
for (auto k = begin; k < end; ++k) r.terms[k].second /= g;
});
}
}
if (!detail::propagate(p.rows, p.vars, budget)) {
p.infeasible = true;
return p;
}
}
p.map = budget.filled<detail::Mapping>(vars_.size());
p.c.reserve(vars_.size());
p.a.reserve(p.rows.size() + p.vars.size());
p.b.reserve(p.rows.size() + p.vars.size());
const double sign = sense_ == Sense::Maximize ? 1 : -1;
p.offset = sign * constant_;
auto new_column = [&](int j, double coef) {
p.map[j].terms.emplace_back(int(p.c.size()), coef);
p.c.push_back(sign * objective_[j] * coef);
};
for (int j = 0; j < int(vars_.size()); ++j) {
budget.checkpoint(1);
const auto &v = p.vars[j];
if (v.lower == v.upper)
p.map[j].offset = v.lower;
else if (std::isfinite(v.lower)) {
p.map[j].offset = v.lower;
new_column(j, 1);
} else if (std::isfinite(v.upper)) {
p.map[j].offset = v.upper;
new_column(j, -1);
} else {
new_column(j, 1);
new_column(j, -1);
}
p.offset += sign * objective_[j] * p.map[j].offset;
}
auto append = [&](Terms a, double b) {
double largest = 0;
budget.chunks(a.size(), [&](std::size_t begin, std::size_t end) {
for (auto k = begin; k < end; ++k)
largest = std::max(largest, std::abs(a[k].second));
});
if (largest == 0) {
if (b < -opt.feasibility_tolerance) p.infeasible = true;
return;
}
if (opt.scaling) {
// 2 の冪で行をスケール。係数比を悪化させず丸めも避ける。
int e = 0;
std::frexp(largest, &e);
double scale = std::ldexp(1.0, -std::max(-500, std::min(500, e)));
budget.chunks(a.size(), [&](std::size_t begin, std::size_t end) {
for (auto k = begin; k < end; ++k) a[k].second *= scale;
});
b *= scale;
}
p.a.push_terms(std::move(a));
p.b.push_back(b);
};
for (const auto &r : p.rows) {
budget.checkpoint(1);
long double rhs = r.rhs;
Terms a;
a.reserve(r.terms.size());
long double maximum = 0;
bool bounded = true;
budget.chunks(r.terms.size(), [&](std::size_t begin, std::size_t end) {
for (auto k = begin; k < end; ++k) {
const auto [j, coef] = r.terms[k];
const auto &v = p.vars[j];
double extreme = coef > 0 ? v.upper : v.lower;
if (std::isfinite(extreme))
maximum += (long double)coef * extreme;
else
bounded = false;
rhs -= (long double)coef * p.map[j].offset;
for (auto [column, scale] : p.map[j].terms)
a.emplace_back(column, coef * scale);
}
});
if (opt.presolve && bounded && maximum <= (long double)r.rhs) continue;
append(std::move(a), double(rhs));
}
// 非負係数の残存行から直接導ける上限は重複して追加しない。
// 伝播で得た上限をすべて行にすると、割当 LP 等で n^2 行増えてしまう。
auto implied_upper = budget.filled<double>(p.c.size(), inf);
if (opt.presolve)
for (int i = 0; i < int(p.a.size()); ++i) {
bool nonnegative = true;
budget.chunks(p.a[i].size(), [&](std::size_t begin, std::size_t end) {
for (auto k = begin; k < end; ++k)
if (p.a[i][k].second < 0) {
nonnegative = false;
break;
}
});
if (!nonnegative) continue;
budget.chunks(p.a[i].size(), [&](std::size_t begin, std::size_t end) {
for (auto k = begin; k < end; ++k) {
auto [j, v] = p.a[i][k];
if (v > 0) implied_upper[j] = std::min(implied_upper[j], p.b[i] / v);
}
});
}
for (int j = 0; j < int(p.vars.size()); ++j) {
budget.checkpoint(1);
const auto &v = p.vars[j];
if (v.lower == v.upper || !std::isfinite(v.lower) || !std::isfinite(v.upper)) continue;
if (implied_upper[p.map[j].terms[0].first] <= v.upper - v.lower) continue;
Terms a;
a.emplace_back(p.map[j].terms[0].first, 1);
append(std::move(a), v.upper - v.lower);
}
return p;
}
// 古い前処理で得た境界と変数変換は、制約を追加しても有効なまま。
// 再度の列消去を避け、追加行だけを保存済みの座標に変換する。
void append_to_root(detail::Prepared &p, detail::Tableau &tab, std::size_t first,
detail::Budget &budget) const {
const auto &opt = budget.opt;
for (std::size_t i = first; i < rows_.size(); ++i) {
budget.checkpoint();
const auto &r = rows_[i];
long long g = 0;
bool integral = opt.presolve;
if (integral)
for (auto [j, coef] : r.terms) {
budget.checkpoint(1);
if (p.vars[j].type == VarType::Continuous || std::abs(coef) > 1e9 ||
coef != std::round(coef)) {
integral = false;
break;
}
g = std::gcd(g, (long long)std::llround(std::abs(coef)));
}
const double divisor = integral && g ? double(g) : 1;
long double rhs =
integral && g ? std::floor(r.rhs / divisor + opt.integrality_tolerance) : r.rhs;
Terms a;
a.reserve(r.terms.size());
budget.chunks(r.terms.size(), [&](std::size_t begin, std::size_t end) {
for (auto k = begin; k < end; ++k) {
const auto [j, raw] = r.terms[k];
const double coef = raw / divisor;
rhs -= (long double)coef * p.map[j].offset;
for (auto [column, multiplier] : p.map[j].terms)
a.emplace_back(column, coef * multiplier);
}
});
double largest = 0;
budget.chunks(a.size(), [&](std::size_t begin, std::size_t end) {
for (auto j = begin; j < end; ++j)
largest = std::max(largest, std::abs(a[j].second));
});
if (!largest) {
if (rhs < -opt.feasibility_tolerance) {
p.infeasible = true;
return;
}
continue;
}
double scale = 1;
if (opt.scaling) {
int e = 0;
std::frexp(largest, &e);
scale = std::ldexp(1.0, -std::max(-500, std::min(500, e)));
}
budget.chunks(a.size(), [&](std::size_t begin, std::size_t end) {
for (auto k = begin; k < end; ++k) a[k].second *= scale;
});
const double b = double(rhs) * scale;
tab.add_constraint(a, b, &budget);
p.a.push_terms(std::move(a));
p.b.push_back(b);
}
}
void save_root(detail::Prepared &p, detail::Tableau &tab, detail::Budget &budget, bool relax,
bool consume) const {
budget.checkpoint();
auto next = std::make_shared<RootCache>();
if (consume) {
// この根で最適性を証明済みなら、直後に返すので辞書の複製は不要。
next->prepared = std::move(p);
next->prepared.rows.clear();
next->tableau = std::move(tab);
} else {
next->prepared = p.clone(budget);
next->tableau = tab.clone(budget);
}
next->rows = rows_.size();
const auto &o = budget.opt;
next->presolve = o.presolve;
next->scaling = o.scaling;
next->dualize = o.dualize;
next->relax = relax;
next->feasibility = o.feasibility_tolerance;
next->integrality = o.integrality_tolerance;
next->pivot = o.pivot_tolerance;
root_cache_ = std::move(next); // 完全なスナップショットができてから差し替える。
}
std::vector<double> recover(const detail::Prepared &p, const std::vector<double> &y) const {
std::vector<double> x(vars_.size());
for (int j = 0; j < int(x.size()); ++j) {
x[j] = p.map[j].offset;
for (auto [k, a] : p.map[j].terms) x[j] += a * y[k];
}
return x;
}
// 元変数を直接復元し、変換後の全列ぶんの一時ベクトルを作らない。
std::vector<double> recover(const detail::Prepared &p, const detail::Tableau &tab) const {
std::vector<double> x(vars_.size());
for (int j = 0; j < int(x.size()); ++j) {
x[j] = p.map[j].offset;
for (auto [k, a] : p.map[j].terms) {
const int pos = tab.location[k];
const double y = pos < 0 ? tab.row(-pos - 1)[tab.n + 1] : 0;
x[j] += a * y;
}
}
return x;
}
double objective_value(const std::vector<double> &x) const {
long double value = constant_;
for (int j = 0; j < int(x.size()); ++j) value += (long double)objective_[j] * x[j];
return double(value);
}
bool
feasible(const std::vector<double> &x, double tol, bool relax, detail::Budget *budget) const {
if (x.size() != vars_.size()) return false;
if (!detail::Budget::scan(budget, x.size(), [&](std::size_t begin, std::size_t end) {
for (auto j = begin; j < end; ++j) {
if (!std::isfinite(x[j]) || x[j] < vars_[j].lower - tol ||
x[j] > vars_[j].upper + tol)
return false;
if (!relax && vars_[j].type != VarType::Continuous &&
std::abs(x[j] - std::round(x[j])) > tol)
return false;
}
return true;
}))
return false;
for (const auto &r : rows_) {
if (budget) budget->checkpoint(1);
long double lhs = 0, magnitude = std::abs(r.rhs);
auto evaluate = [&](std::size_t begin, std::size_t end) {
for (auto k = begin; k < end; ++k) {
const auto [j, a] = r.terms[k];
lhs += (long double)a * x[j];
magnitude += std::abs((long double)a * x[j]);
}
};
if (budget)
budget->chunks(r.terms.size(), evaluate);
else
evaluate(0, r.terms.size());
if (lhs > (long double)r.rhs + tol + 1e-12L * magnitude) return false;
}
return true;
}
public:
int add_variable(double lower = 0, double upper = inf, VarType type = VarType::Continuous) {
if (std::isnan(lower) || std::isnan(upper) || lower == inf || upper == -inf)
throw std::invalid_argument("invalid variable bound");
if (type == VarType::Binary) {
lower = std::max(lower, 0.0);
upper = std::min(upper, 1.0);
}
root_cache_.reset();
vars_.push_back({lower, upper, type});
objective_.push_back(0);
return int(vars_.size()) - 1;
}
int add_integer(double lower = 0, double upper = inf) {
return add_variable(lower, upper, VarType::Integer);
}
int add_binary() { return add_variable(0, 1, VarType::Binary); }
int num_variables() const { return int(vars_.size()); }
void clear_root_cache() const { root_cache_.reset(); }
void add_constraint(Terms terms, Relation rel, double rhs) {
if (!std::isfinite(rhs)) throw std::invalid_argument("non-finite rhs");
terms = normalize(std::move(terms));
if (rel == Relation::LessEqual) {
rows_.push_back({std::move(terms), rhs});
return;
}
if (rel != Relation::GreaterEqual) rows_.push_back({terms, rhs});
if (rel != Relation::LessEqual) {
for (auto &t : terms) t.second = -t.second;
rows_.push_back({std::move(terms), -rhs});
}
}
void add_le(Terms a, double b) { add_constraint(std::move(a), Relation::LessEqual, b); }
void add_ge(Terms a, double b) { add_constraint(std::move(a), Relation::GreaterEqual, b); }
void add_eq(Terms a, double b) { add_constraint(std::move(a), Relation::Equal, b); }
void set_objective(Terms terms, Sense sense = Sense::Minimize, double constant = 0) {
if (!std::isfinite(constant)) throw std::invalid_argument("non-finite objective constant");
terms = normalize(std::move(terms));
root_cache_.reset();
std::fill(objective_.begin(), objective_.end(), 0);
for (auto [j, a] : terms) objective_[j] = a;
sense_ = sense;
constant_ = constant;
}
void minimize(Terms terms, double constant = 0) {
set_objective(std::move(terms), Sense::Minimize, constant);
}
void maximize(Terms terms, double constant = 0) {
set_objective(std::move(terms), Sense::Maximize, constant);
}
bool is_feasible(const std::vector<double> &x, double tol = 1e-7, bool relax = false) const {
return feasible(x, tol, relax, nullptr);
}
// relax=true なら整数変数も連続変数として解く。
Result solve(const Options &opt = {}, bool relax = false) const {
if (!(opt.time_limit >= 0) || !(opt.feasibility_tolerance > 0) ||
!(opt.integrality_tolerance > 0) || !(opt.pivot_tolerance > 0) ||
!(opt.absolute_gap >= 0) || !(opt.relative_gap >= 0) ||
!std::isfinite(opt.feasibility_tolerance) ||
!std::isfinite(opt.integrality_tolerance) || !std::isfinite(opt.pivot_tolerance) ||
!std::isfinite(opt.absolute_gap) || !std::isfinite(opt.relative_gap) ||
opt.strong_branching_candidates < 0 || opt.strong_branching_iterations < 0 ||
opt.strong_branching_depth < 0 || opt.root_cut_rounds < -1 || opt.cuts_per_round < 0)
throw std::invalid_argument("invalid solver options");
Result ret;
detail::Budget budget{opt, ret.stats};
ret.stats.original_variables = int(vars_.size());
ret.stats.original_constraints = int(rows_.size());
const double sign = sense_ == Sense::Maximize ? 1 : -1;
double incumbent = -inf;
double global_bound = inf;
auto finish = [&](Status status) {
if (status == Status::Infeasible && ret.has_solution) status = Status::NumericalError;
// カットで条件が悪くなった場合は元のモデルへ戻る。暫定解と総予算を継承。
if (status == Status::NumericalError && ret.stats.cuts) {
Options retry = opt;
retry.root_cut_rounds = 0;
retry.time_limit = std::max(0.0, opt.time_limit - budget.elapsed());
retry.iteration_limit =
opt.iteration_limit - std::min(opt.iteration_limit, ret.stats.pivots);
retry.node_limit = opt.node_limit - std::min(opt.node_limit, ret.stats.nodes);
if (ret.has_solution) retry.initial_solution = ret.x;
Result fresh = solve(retry, relax);
fresh.stats.nodes += ret.stats.nodes;
fresh.stats.pivots += ret.stats.pivots;
fresh.stats.strong_branching_pivots += ret.stats.strong_branching_pivots;
fresh.stats.root_pivots += ret.stats.root_pivots;
fresh.stats.reused_root = fresh.stats.reused_root || ret.stats.reused_root;
fresh.stats.heuristic_solutions += ret.stats.heuristic_solutions;
fresh.stats.cuts += ret.stats.cuts;
fresh.stats.reduced_cost_fixings += ret.stats.reduced_cost_fixings;
fresh.stats.numerical_restarts += ret.stats.numerical_restarts + 1;
fresh.stats.seconds = budget.elapsed();
return fresh;
}
ret.status = status;
if (ret.has_solution) global_bound = std::max(global_bound, incumbent);
ret.best_bound = sign * global_bound;
ret.stats.seconds = budget.elapsed();
return ret;
};
auto remember = [&](const std::vector<double> &x, bool heuristic) {
double val = objective_value(x) * sign;
if (!std::isfinite(val) || val <= incumbent) return false;
incumbent = val;
ret.has_solution = true;
ret.x = x;
ret.objective = sign * val;
if (heuristic) ++ret.stats.heuristic_solutions;
return true;
};
auto accept = [&](const std::vector<double> &x, bool heuristic, bool timed_check = true) {
if (!feasible(x, opt.feasibility_tolerance, relax, timed_check ? &budget : nullptr))
return false;
return remember(x, heuristic);
};
try {
if (!opt.initial_solution.empty()) {
if (!is_feasible(opt.initial_solution, opt.feasibility_tolerance, relax))
throw std::invalid_argument("initial_solution is not feasible");
remember(opt.initial_solution, true);
}
if (budget.stop() || opt.node_limit == 0) return finish(Status::Limit);
auto cached = opt.reuse_root && root_cache_ && root_cache_->matches(opt, relax)
? root_cache_
: nullptr;
bool reuse = bool(cached);
detail::Prepared p;
detail::Tableau root;
std::vector<int> integers;
std::vector<double> root_x;
for (;;) {
if (budget.stop() || ret.stats.nodes >= opt.node_limit)
return finish(Status::Limit);
if (reuse) {
ret.stats.reused_root = true;
global_bound = cached->tableau.value() + cached->prepared.offset;
p = cached->prepared.clone(budget);
root = cached->tableau.clone(budget);
append_to_root(p, root, cached->rows, budget);
} else
p = prepare(budget, relax);
budget.checkpoint();
if (p.infeasible) {
global_bound = -inf;
return finish(Status::Infeasible);
}
if (opt.rounding_heuristic) {
// 根 LP が時間切れでも使える簡単な暫定解を先に検査する。
for (bool upper : {false, true}) {
auto x = budget.filled<double>(p.vars.size());
budget.chunks(p.vars.size(), [&](std::size_t begin, std::size_t end) {
for (auto j = begin; j < end; ++j) {
const auto &v = p.vars[j];
double target = upper ? v.upper : v.lower;
x[j] = std::isfinite(target)
? target
: std::max(v.lower, std::min(v.upper, 0.0));
}
});
accept(x, true);
}
}
ret.stats.presolved_variables = int(p.c.size());
ret.stats.presolved_constraints = int(p.b.size());
integers.clear();
for (int j = 0; j < int(vars_.size()); ++j)
if (!relax && vars_[j].type != VarType::Continuous && !p.map[j].terms.empty())
integers.push_back(j);
// 多数の非負変数を持つ covering/assignment LP は、双対のほうが小さい
// 基底で Phase I も不要になることが多い。MIP は分枝用の主辞書を保持する。
if (!reuse && opt.dualize && integers.empty() && !p.b.empty() &&
p.c.size() > 2 * p.b.size() &&
std::all_of(p.c.begin(), p.c.end(), [](double c) { return c <= 0; })) {
std::vector<std::vector<double>> a;
a.reserve(p.c.size());
for (std::size_t j = 0; j < p.c.size(); ++j)
a.push_back(budget.filled<double>(p.b.size()));
std::vector<double> b = p.c, c = p.b;
for (double &v : b) v = -v;
for (double &v : c) v = -v;
for (int i = 0; i < int(p.b.size()); ++i)
budget.chunks(p.a[i].size(), [&](std::size_t begin, std::size_t end) {
for (auto k = begin; k < end; ++k) {
auto [j, v] = p.a[i][k];
a[j][i] = -v;
}
});
detail::Tableau dual(a, b, c, &budget);
++ret.stats.nodes;
const auto before = ret.stats.pivots;
const Status s = dual.solve(budget);
ret.stats.root_pivots += ret.stats.pivots - before;
if (s == Status::Limit) return finish(s);
if (s == Status::Optimal) {
// 双対の各制約の slack に対する reduced cost が元 LP の変数値。
std::vector<double> y(p.c.size());
for (int j = 0; j < int(y.size()); ++j) {
const int pos = dual.location[dual.n + j];
if (pos >= 0) y[j] = dual.row(dual.m)[pos] / dual.objective_scale;
}
auto x = recover(p, y);
const double obj = sign * objective_value(x);
const double bound = p.offset - dual.value();
if (std::isfinite(obj) && std::isfinite(bound) &&
is_feasible(x, opt.feasibility_tolerance, relax) &&
std::abs(obj - bound) <=
opt.feasibility_tolerance * (1 + std::abs(obj))) {
remember(x, false);
global_bound = incumbent;
ret.stats.solved_dual = true;
return finish(Status::Optimal);
}
}
// 双対の状態だけでは曖昧な場合や復元誤差がある場合は、主問題で判定。
}
if (budget.stop() || ret.stats.nodes >= opt.node_limit)
return finish(Status::Limit);
if (!reuse) root = detail::Tableau(p.a, p.b, p.c, &budget);
++ret.stats.nodes;
const auto before = ret.stats.pivots;
// 長い再最適化は一度だけ作り直す。総時間・pivot・node の予算は継続。
const int soft_limit =
int(std::min<long long>(std::numeric_limits<int>::max(),
std::max(256LL, 4LL * (root.n + (long long)root.m))));
Status root_status = reuse ? root.dual(budget, soft_limit) : root.solve(budget);
ret.stats.root_pivots += ret.stats.pivots - before;
if (root_status == Status::Optimal) {
root_x = recover(p, root);
const double obj = sign * objective_value(root_x);
const double bound = root.value() + p.offset;
bool valid =
std::isfinite(obj) && std::isfinite(bound) &&
feasible(root_x, opt.feasibility_tolerance, true, &budget) &&
std::abs(obj - bound) <= opt.feasibility_tolerance * (1 + std::abs(obj));
for (int j = 0; j <= root.n && valid; ++j)
if (root.nonbasic[j] != -1)
valid = std::isfinite(root.row(root.m)[j]) &&
root.row(root.m)[j] >= -10 * opt.pivot_tolerance;
if (!valid) root_status = Status::NumericalError;
}
if (reuse && root_status != Status::Optimal &&
(root_status != Status::Limit || !budget.stop())) {
// 不変の根以外は保存しない。再利用中の不整合・実行不能も新規求解で確認。
root_cache_.reset();
cached.reset();
reuse = false;
global_bound = inf;
++ret.stats.numerical_restarts;
continue;
}
if (root_status != Status::Optimal) {
// Phase II の打切りなら現在の基底も実行可能候補。整数性も検査する。
if (root_status == Status::Limit && root.primal_feasible)
accept(recover(p, root), false, false);
if (root_status == Status::Infeasible) global_bound = -inf;
if (root_status == Status::Unbounded &&
!is_feasible(recover(p, root), opt.feasibility_tolerance, true))
root_status = Status::NumericalError;
if (root_status == Status::Unbounded && !integers.empty())
root_status = Status::RelaxationUnbounded;
return finish(root_status);
}
break;
}
global_bound = root.value() + p.offset;
// 保存の途中で時間切れになっても、既に確認できた整数解は失わない。
bool integral_root = true;
auto rounded = root_x;
for (int j : integers) {
if (std::abs(rounded[j] - std::round(rounded[j])) > opt.integrality_tolerance)
integral_root = false;
rounded[j] = std::round(rounded[j]);
}
bool root_proven = false;
if (integral_root && accept(rounded, false))
root_proven = std::abs(incumbent - global_bound) <=
opt.feasibility_tolerance * (1 + std::abs(incumbent));
// 同じ値の initial_solution が既にあっても根で最適性を証明できる。
if (integral_root && ret.has_solution &&
std::abs(incumbent - global_bound) <=
opt.feasibility_tolerance * (1 + std::abs(incumbent)))
root_proven = true;
if (opt.reuse_root && (!reuse || cached->rows != rows_.size())) {
try {
save_root(p, root, budget, relax, root_proven);
} catch (const detail::TimeLimitReached &) {
if (!root_proven) throw;
}
}
if (root_proven) {
global_bound = incumbent;
return finish(Status::Optimal);
}
if (integers.empty()) return finish(Status::NumericalError);
std::vector<double> rounding_buffer;
auto round_candidate = [&](const std::vector<double> &x) {
for (int mode = 0; mode < 3; ++mode) {
rounding_buffer.assign(x.begin(), x.end());
auto &y = rounding_buffer;
for (int j : integers)
y[j] = mode == 0 ? std::round(y[j])
: mode == 1 ? std::floor(y[j])
: std::ceil(y[j]);
accept(y, true);
}
};
if (opt.rounding_heuristic) round_candidate(root_x);
if (opt.rounding_heuristic && ret.has_solution) {
// 実行可能な丸め解から、目的を改善する整数の1座標移動を試す。
// 問題の形は仮定せず、各候補を元の全制約で検査する。
auto order = budget.copy(integers);
std::sort(order.begin(), order.end(), [&](int a, int b) {
budget.checkpoint(1);
return std::abs(objective_[a]) > std::abs(objective_[b]);
});
auto x = budget.copy(ret.x);
for (int j : order) {
budget.checkpoint(1);
const double gain = sign * objective_[j];
if (gain == 0) continue;
const double old = x[j], next = old + (gain > 0 ? 1 : -1);
if (next < p.vars[j].lower || next > p.vars[j].upper) continue;
x[j] = next;
if (feasible(x, opt.feasibility_tolerance, false, &budget))
remember(x, true);
else
x[j] = old;
}
}
std::vector<bool> integer_column(p.c.size(), false);
for (int j : integers)
if (p.map[j].terms.size() == 1 && p.map[j].offset == std::round(p.map[j].offset))
integer_column[p.map[j].terms[0].first] = true;
int cut_rounds = opt.root_cut_rounds;
if (cut_rounds < 0) {
const auto x = recover(p, root);
int fractional = 0;
for (int j : integers)
fractional += std::abs(x[j] - std::round(x[j])) > opt.integrality_tolerance;
cut_rounds = fractional == 1 ? 3 : 0;
}
for (int round = 0; round < cut_rounds; ++round) {
if (budget.stop()) return finish(Status::Limit);
const int added = detail::gomory_cuts(root, p, integer_column, budget);
if (!added) break;
ret.stats.cuts += added;
Status s = root.dual(budget);
if (s == Status::Infeasible) {
global_bound = -inf;
return finish(s);
}
if (s != Status::Optimal) return finish(s);
global_bound = root.value() + p.offset;
}
// 整数目的関数なら LP 上界を格子に切り下げ、無駄な探索を省く。
double lattice_base = sign * constant_, lattice_step = 0;
if (opt.integer_objective) {
long long g = 0;
bool valid = true;
for (int j = 0; j < int(vars_.size()); ++j) {
if (p.map[j].terms.empty()) {
lattice_base += sign * objective_[j] * p.map[j].offset;
continue;
}
double a = objective_[j];
if (a == 0) continue;
if (vars_[j].type == VarType::Continuous || a != std::round(a) ||
std::abs(a) > 1e9) {
valid = false;
break;
}
g = std::gcd(g, (long long)std::llround(std::abs(a)));
}
if (valid) lattice_step = double(g);
}
auto tighten = [&](double bound) {
if (lattice_step && std::isfinite(bound)) {
const double slack =
opt.integrality_tolerance +
1e-12 * (std::abs(bound) + std::abs(lattice_base)) / lattice_step;
return lattice_base +
lattice_step * std::floor((bound - lattice_base) / lattice_step + slack);
}
return bound;
};
auto prunable = [&](double b) {
return ret.has_solution && tighten(b) <= incumbent + opt.absolute_gap +
opt.relative_gap * std::abs(incumbent);
};
if (opt.rounding_heuristic) round_candidate(recover(p, root));
if (prunable(global_bound)) {
global_bound = tighten(global_bound);
return finish(Status::Optimal);
}
auto improvement_target = [&] {
if (!lattice_step) return incumbent;
return lattice_base +
lattice_step * (std::floor((incumbent - lattice_base) / lattice_step) + 1);
};
std::vector<double> step;
if (opt.reduced_cost_fixing) {
step = budget.filled<double>(root.next_id);
for (int j = 0; j < root.original_n; ++j)
if (integer_column[j]) step[j] = 1;
// 整数変数の上下限に対応する slack も、スケーリング後の刻みを持つ。
for (std::size_t i = 0; i < p.a.size(); ++i) {
budget.checkpoint(1);
if (p.a[i].size() != 1 || !integer_column[p.a[i][0].first]) continue;
const double scale = std::abs(p.a[i][0].second);
const double rhs = p.b[i] / scale;
if (std::isfinite(rhs) && rhs == std::round(rhs))
step[root.original_n + i] = scale;
}
if (ret.has_solution)
ret.stats.reduced_cost_fixings +=
root.fix_nonbasic(step, global_bound - improvement_target(), budget);
}
struct Node {
detail::Tableau tab;
std::vector<detail::Row> branches;
double bound;
int depth;
bool solved;
int branch_var = -1;
bool up = false;
double distance = 0, parent_bound = 0;
};
std::vector<Node> pending;
pending.push_back({std::move(root), {}, global_bound, 0, true});
struct Pseudo {
double down = 0, up = 0;
int nd = 0, nu = 0;
};
std::vector<Pseudo> pc(vars_.size());
auto add_branch = [&](detail::Tableau &tab, int j, bool up, double bound) {
budget.checkpoint();
const double s = up ? -1 : 1;
tab.add_constraint(p.map[j].terms, s * (bound - p.map[j].offset), nullptr, s);
};
auto update_pc = [&](int j, bool up, double distance, double gain) {
if (!std::isfinite(gain)) return;
auto &v = pc[j];
(up ? v.up : v.down) += std::max(0.0, gain) / std::max(distance, 1e-6);
++(up ? v.nu : v.nd);
};
std::vector<std::pair<double, int>> candidates;
while (!pending.empty()) {
global_bound = incumbent;
for (const auto &node : pending)
global_bound = std::max(global_bound, tighten(node.bound));
if (prunable(global_bound)) return finish(Status::Optimal);
if (budget.stop()) return finish(Status::Limit);
Node node = std::move(pending.back());
pending.pop_back();
if (prunable(node.bound)) continue;
if (node.depth) {
if (ret.stats.nodes >= opt.node_limit) return finish(Status::Limit);
++ret.stats.nodes;
}
if (!node.solved) {
Status s;
if (opt.warm_start)
s = node.tab.dual(budget);
else {
node.tab = detail::Tableau(p.a, p.b, p.c, &budget);
// コールド比較用: 分枝行を元の非基底辞書に追加し Phase I から解く。
for (const auto &b : node.branches) {
int j = b.terms[0].first;
bool up = b.terms[0].second < 0;
add_branch(node.tab, j, up, up ? -b.rhs : b.rhs);
node.tab.write_row(node.tab.m - 1)[node.tab.n] = -1;
}
s = node.tab.solve(budget);
}
if (s == Status::Infeasible) continue;
if (s != Status::Optimal) return finish(s);
node.bound = node.tab.value() + p.offset;
if (node.branch_var >= 0)
update_pc(node.branch_var, node.up, node.distance,
node.parent_bound - node.bound);
if (prunable(node.bound)) continue;
}
auto x = recover(p, node.tab);
if (!is_feasible(x, opt.feasibility_tolerance, true))
return finish(Status::NumericalError);
if (opt.rounding_heuristic) round_candidate(x);
if (prunable(node.bound)) continue;
if (opt.reduced_cost_fixing && ret.has_solution && node.depth) {
step.resize(node.tab.next_id, 1);
ret.stats.reduced_cost_fixings +=
node.tab.fix_nonbasic(step, node.bound - improvement_target(), budget,
std::max(4, node.tab.n / 8));
}
candidates.clear();
for (int j : integers) {
const double f = x[j] - std::floor(x[j]);
if (std::min(f, 1 - f) <= opt.integrality_tolerance) continue;
const auto &v = pc[j];
double d = f, u = 1 - f;
if (opt.pseudocost) {
d *= v.nd ? v.down / v.nd : 1;
u *= v.nu ? v.up / v.nu : 1;
}
candidates.emplace_back(std::min(d, u) + 0.1 * std::max(d, u), j);
}
if (candidates.empty()) {
for (int j : integers) x[j] = std::round(x[j]);
if (!is_feasible(x, opt.feasibility_tolerance, false))
return finish(Status::NumericalError);
remember(x, false);
continue;
}
std::sort(
candidates.begin(), candidates.end(), std::greater<std::pair<double, int>>());
int branch = candidates[0].second;
double best_score = -1;
detail::Tableau best_down, best_up;
Status down_status = Status::Limit, up_status = Status::Limit;
bool probed = false;
const int probes =
node.depth < opt.strong_branching_depth
? std::min(int(candidates.size()), opt.strong_branching_candidates)
: 0;
for (int k = 0; k < probes; ++k) {
int j = candidates[k].second;
const double f = x[j] - std::floor(x[j]);
budget.checkpoint();
detail::Tableau down = node.tab, up = node.tab;
add_branch(down, j, false, std::floor(x[j]));
add_branch(up, j, true, std::ceil(x[j]));
auto before = ret.stats.pivots;
auto ds = down.dual(budget, opt.strong_branching_iterations);
auto us = up.dual(budget, opt.strong_branching_iterations);
ret.stats.strong_branching_pivots += ret.stats.pivots - before;
if (budget.stop()) return finish(Status::Limit);
if (ds == Status::NumericalError || us == Status::NumericalError)
return finish(Status::NumericalError);
double dg = ds == Status::Infeasible
? inf
: std::max(0.0, node.bound - (down.value() + p.offset));
double ug = us == Status::Infeasible
? inf
: std::max(0.0, node.bound - (up.value() + p.offset));
if (ds == Status::Optimal) update_pc(j, false, f, dg);
if (us == Status::Optimal) update_pc(j, true, 1 - f, ug);
if (opt.rounding_heuristic) {
if (ds == Status::Optimal) round_candidate(recover(p, down));
if (us == Status::Optimal) round_candidate(recover(p, up));
}
const double score = std::min(dg, ug) + 0.1 * std::max(dg, ug);
if (!probed || score > best_score) {
probed = true;
branch = j;
best_score = score;
best_down = std::move(down);
best_up = std::move(up);
down_status = ds;
up_status = us;
}
if (ds == Status::Infeasible || us == Status::Infeasible) break;
}
if (!probed) {
budget.checkpoint();
best_down = node.tab;
best_up = std::move(node.tab);
add_branch(best_down, branch, false, std::floor(x[branch]));
add_branch(best_up, branch, true, std::ceil(x[branch]));
}
auto push = [&](detail::Tableau tab, bool up, Status s) {
if (s == Status::Infeasible) return;
double b = std::min(node.bound, tab.value() + p.offset);
if (prunable(b)) return;
std::vector<detail::Row> branches;
if (!opt.warm_start) {
branches = node.branches;
branches.push_back({{{branch, up ? -1.0 : 1.0}},
up ? -std::ceil(x[branch]) : std::floor(x[branch])});
}
const double distance =
up ? std::ceil(x[branch]) - x[branch] : x[branch] - std::floor(x[branch]);
pending.push_back({std::move(tab), std::move(branches), b, node.depth + 1,
s == Status::Optimal, branch, up, distance, node.bound});
};
// 有望な側を先に深さ優先探索し、暫定解を早めに得る。
bool up_first = best_up.value() > best_down.value();
if (up_first) {
push(std::move(best_down), false, down_status);
push(std::move(best_up), true, up_status);
} else {
push(std::move(best_up), true, up_status);
push(std::move(best_down), false, down_status);
}
}
global_bound = incumbent;
return finish(ret.has_solution ? Status::Optimal : Status::Infeasible);
} catch (const detail::TimeLimitReached &) { return finish(Status::Limit); }
}
Result solve_relaxation(const Options &opt = {}) const { return solve(opt, true); }
};
// 同じ変数で制約追加・目的関数変更を繰り返す LP 向け。
// maximize c*x, A*x<=b, x>=0。整数制約・変数変換・presolve は Model を使う。
// 制約追加だけなら双対法、目的関数変更だけなら主単体法で基底を再利用する。
class LPSolver {
enum class State { Cold, Optimal, Primal, Dual };
State state_ = State::Cold;
std::vector<std::vector<double>> a_;
std::vector<double> b_, c_;
detail::Tableau tab_;
void check_coefficients(const std::vector<double> &row) const {
if (row.size() != c_.size()) throw std::invalid_argument("LP dimension mismatch");
for (double v : row)
if (!std::isfinite(v)) throw std::invalid_argument("non-finite LP coefficient");
}
public:
LPSolver(std::vector<std::vector<double>> a, std::vector<double> b, std::vector<double> c)
: a_(std::move(a)), b_(std::move(b)), c_(std::move(c)) {
if (a_.size() != b_.size()) throw std::invalid_argument("LP A/b dimension mismatch");
check_coefficients(c_);
for (const auto &row : a_) check_coefficients(row);
for (double v : b_)
if (!std::isfinite(v)) throw std::invalid_argument("non-finite LP rhs");
}
void add_le(std::vector<double> a, double b) {
check_coefficients(a);
if (!std::isfinite(b)) throw std::invalid_argument("non-finite LP rhs");
if (state_ == State::Optimal || state_ == State::Dual) {
Terms terms;
for (int j = 0; j < int(a.size()); ++j)
if (a[j]) terms.emplace_back(j, a[j]);
tab_.add_constraint(terms, b);
state_ = State::Dual;
} else
state_ = State::Cold;
a_.push_back(std::move(a));
b_.push_back(b);
}
void set_objective(std::vector<double> c) {
check_coefficients(c);
c_ = std::move(c);
if (state_ == State::Optimal || state_ == State::Primal) {
tab_.scale_objective(c_);
double *z = tab_.write_row(tab_.m);
std::fill(z, z + tab_.stride, 0);
for (int j = 0; j <= tab_.n; ++j) {
int id = tab_.nonbasic[j];
if (id >= 0 && id < tab_.n) z[j] = -c_[id] * tab_.objective_scale;
}
for (int i = 0; i < tab_.m; ++i) {
int id = tab_.basic[i];
if (id < 0 || id >= tab_.n || c_[id] == 0) continue;
const double *row = tab_.row(i);
for (int j = 0; j < tab_.stride; ++j)
z[j] += c_[id] * tab_.objective_scale * row[j];
}
state_ = State::Primal;
} else
state_ = State::Cold;
}
Result solve(const Options &opt = {}) {
if (!(opt.time_limit >= 0) || !(opt.feasibility_tolerance > 0) ||
!(opt.pivot_tolerance > 0) || !std::isfinite(opt.feasibility_tolerance) ||
!std::isfinite(opt.pivot_tolerance))
throw std::invalid_argument("invalid LP options");
Result r;
r.best_bound = inf;
r.stats.original_variables = r.stats.presolved_variables = int(c_.size());
r.stats.original_constraints = r.stats.presolved_constraints = int(b_.size());
detail::Budget budget{opt, r.stats};
try {
if (budget.stop() || opt.node_limit == 0)
r.status = Status::Limit;
else {
++r.stats.nodes;
if (state_ == State::Cold) {
tab_ = detail::Tableau(a_, b_, c_, &budget);
state_ = State::Primal;
r.status = tab_.solve(budget);
} else if (state_ == State::Dual)
r.status = tab_.dual(budget);
else if (state_ == State::Primal)
r.status = tab_.primal(tab_.m, budget);
else
r.status = Status::Optimal;
}
} catch (const detail::TimeLimitReached &) { r.status = Status::Limit; }
if (r.status == Status::Optimal || (r.status == Status::Limit && tab_.primal_feasible &&
state_ != State::Cold && state_ != State::Dual)) {
r.x = tab_.solution();
bool valid = true;
long double objective = 0;
for (int j = 0; j < int(c_.size()); ++j) {
if (!std::isfinite(r.x[j]) || r.x[j] < -opt.feasibility_tolerance) valid = false;
objective += (long double)c_[j] * r.x[j];
}
for (int i = 0; i < int(b_.size()); ++i) {
long double lhs = 0, magnitude = std::abs(b_[i]);
for (int j = 0; j < int(c_.size()); ++j) {
lhs += (long double)a_[i][j] * r.x[j];
magnitude += std::abs((long double)a_[i][j] * r.x[j]);
}
if (lhs > (long double)b_[i] + opt.feasibility_tolerance + 1e-12L * magnitude)
valid = false;
}
if (std::abs(double(objective) - tab_.value()) >
opt.feasibility_tolerance * (1 + std::abs(double(objective))))
valid = false;
if (valid && std::isfinite(double(objective))) {
r.has_solution = true;
r.objective = double(objective);
if (r.status == Status::Optimal) r.best_bound = r.objective;
} else {
if (r.status == Status::Optimal) r.status = Status::NumericalError;
r.x.clear();
}
}
if (r.status == Status::Infeasible) r.best_bound = -inf;
state_ = r.status == Status::Optimal ? State::Optimal : State::Cold;
r.stats.seconds = budget.elapsed();
return r;
}
};
// 行列形式の簡単な入口。maximize c*x subject to A*x <= b, x >= 0。
inline Result linprog(const std::vector<std::vector<double>> &a, const std::vector<double> &b,
const std::vector<double> &c, const Options &options = {}) {
if (a.size() != b.size()) throw std::invalid_argument("A/b dimension mismatch");
Model model;
for (std::size_t j = 0; j < c.size(); ++j) model.add_variable();
Terms cost;
for (int j = 0; j < int(c.size()); ++j) cost.emplace_back(j, c[j]);
model.maximize(std::move(cost));
for (int i = 0; i < int(a.size()); ++i) {
if (a[i].size() != c.size()) throw std::invalid_argument("A/c dimension mismatch");
Terms row;
for (int j = 0; j < int(c.size()); ++j)
if (a[i][j] != 0) row.emplace_back(j, a[i][j]);
model.add_le(std::move(row), b[i]);
}
return model.solve(options);
}
} // namespace cpopt
#endif // mip_HPP
// https://hitonanode.github.io/cplib-cpp/utilities/transpose.hpp
#include <string>
#include <vector>
template <typename T> std::vector<std::vector<T>> transpose(const std::vector<std::vector<T>> &v) {
const int H = v.size(), W = v[0].size();
std::vector<std::vector<T>> ret(W, std::vector<T>(H));
for (int i = 0; i < H; i++) {
for (int j = 0; j < W; j++) ret[j][i] = v[i][j];
}
return ret;
}
std::vector<std::string> transpose(const std::vector<std::string> &v) {
const int H = v.size(), W = v[0].size();
std::vector<std::string> ret(W, std::string(H, '\0'));
for (int i = 0; i < H; i++) {
for (int j = 0; j < W; j++) ret[j][i] = v[i][j];
}
return ret;
}
int main() {
int T;
cin >> T;
while (T--) {
int N, M;
cin >> N >> M;
vector A(N, vector<double>(M));
cin >> A;
// A = transpose(A);
dbg(A);
for (auto &a : A) {
// for (auto &x : a) x = -x;
a.push_back(-1);
a.push_back(1);
}
A.push_back({});
REP(_, A.front().size() - 2) A.back().push_back(1);
A.back().push_back(0);
A.back().push_back(0);
A.push_back({});
REP(_, A.front().size() - 2) A.back().push_back(-1);
A.back().push_back(0);
A.back().push_back(0);
vector<double> B(A.size() - 2, 0);
B.push_back(1);
B.push_back(-1);
vector<double> C(A.front().size());
C.at(C.size() - 2) = -1;
C.back() = 1;
dbg(A);
dbg(B);
dbg(C);
auto ret = cpopt::linprog(A, B, C);
cout << -ret.objective << '\n';
}
}