simplex.hpp¶
求线性规划 \(\max c\mathbin{\cdot}x\)，约束 \(Ax\le b\)、\(x\ge0\)，并区分最优、无解和无界。
Complexity: Time: Exponential worst case; typically polynomial pivot work per iteration. Space: O(rows * columns).
Implementation¶
当前头文件，省略 include guard；依赖见 #include。
/// @complexity Time: Exponential worst case; typically polynomial pivot work per iteration.
/// Space: O(rows * columns).
#include <algorithm>
#include <cassert>
#include <cmath>
#include <limits>
#include <utility>
#include <vector>
namespace noya {
enum class linear_program_status { optimal, infeasible, unbounded };
struct linear_program_result {
linear_program_status sta = linear_program_status::infeasible;
long double obj = 0;
std::vector<long double> sol;
};
namespace simplex_internal {
struct tableau {
static constexpr long double eps = 1e-12L;
// m constraints, n variables.
int m;
int n;
std::vector<int> bas;
std::vector<int> nb;
std::vector<std::vector<long double>> val;
tableau(const std::vector<std::vector<long double>> &a,
const std::vector<long double> &b, const std::vector<long double> &c)
: m(int(b.size())), n(int(c.size())), bas(m), nb(n + 1),
val(m + 2, std::vector<long double>(n + 2)) {
for (int row = 0; row < m; row++) {
for (int col = 0; col < n; col++) {
val[row][col] = a[row][col];
}
bas[row] = n + row;
val[row][n] = -1;
val[row][n + 1] = b[row];
}
for (int col = 0; col < n; col++) {
nb[col] = col;
val[m][col] = -c[col];
}
nb[n] = -1;
val[m + 1][n] = 1;
}
void pivot(int lr, int ec) {
long double inv = 1 / val[lr][ec];
for (int row = 0; row < m + 2; row++) {
if (row == lr) {
continue;
}
for (int col = 0; col < n + 2; col++) {
if (col != ec) {
val[row][col] -= val[lr][col] * val[row][ec] * inv;
}
}
}
for (int col = 0; col < n + 2; col++) {
if (col != ec) {
val[lr][col] *= inv;
}
}
for (int row = 0; row < m + 2; row++) {
if (row != lr) {
val[row][ec] *= -inv;
}
}
val[lr][ec] = inv;
std::swap(bas[lr], nb[ec]);
}
bool optimize(int ph) {
int oi = ph == 1 ? m + 1 : m;
while (true) {
int ent = -1;
for (int col = 0; col <= n; col++) {
if (ph == 2 && nb[col] == -1) {
continue;
}
if (ent == -1 || val[oi][col] < val[oi][ent] - eps ||
(std::abs(val[oi][col] - val[oi][ent]) <= eps &&
nb[col] < nb[ent])) {
ent = col;
}
}
if (ent == -1) {
return true;
}
if (val[oi][ent] >= -eps) {
return true;
}
int lea = -1;
for (int row = 0; row < m; row++) {
if (val[row][ent] <= eps) {
continue;
}
if (lea == -1) {
lea = row;
continue;
}
long double lhs = val[row][n + 1] / val[row][ent];
long double rhs = val[lea][n + 1] / val[lea][ent];
if (lhs < rhs - eps ||
(std::abs(lhs - rhs) <= eps && bas[row] < bas[lea])) {
lea = row;
}
}
if (lea == -1) {
return false;
}
pivot(lea, ent);
}
}
linear_program_result solve() {
int mi = 0;
for (int row = 1; row < m; row++) {
if (val[row][n + 1] < val[mi][n + 1]) {
mi = row;
}
}
if (m > 0 && val[mi][n + 1] < -eps) {
pivot(mi, n);
if (!optimize(1) || val[m + 1][n + 1] < -eps) {
return {linear_program_status::infeasible, 0, {}};
}
if (std::abs(val[m + 1][n + 1]) > eps) {
return {linear_program_status::infeasible, 0, {}};
}
for (int row = 0; row < m; row++) {
if (bas[row] != -1) {
continue;
}
int ent = 0;
for (int col = 1; col <= n; col++) {
if (std::abs(val[row][col]) > std::abs(val[row][ent]) + eps ||
(std::abs(std::abs(val[row][col]) - std::abs(val[row][ent])) <=
eps &&
nb[col] < nb[ent])) {
ent = col;
}
}
if (std::abs(val[row][ent]) > eps) {
pivot(row, ent);
}
}
}
if (!optimize(2)) {
return {linear_program_status::unbounded,
std::numeric_limits<long double>::infinity(),
{}};
}
std::vector<long double> sol(n);
for (int row = 0; row < m; row++) {
if (bas[row] < n) {
sol[bas[row]] = val[row][n + 1];
}
}
return {linear_program_status::optimal, val[m][n + 1], std::move(sol)};
}
};
} // namespace simplex_internal
/// @brief Maximize c*x subject to A*x <= b and x >= 0 with a two-phase
/// simplex tableau; reports optimal, infeasible, or unbounded.
inline linear_program_result
simplex(const std::vector<std::vector<long double>> &a,
const std::vector<long double> &b, const std::vector<long double> &c) {
assert(a.size() == b.size());
for (const auto &row : a) {
assert(row.size() == c.size());
}
return simplex_internal::tableau(a, b, c).solve();
}
} // namespace noya
#ifndef NOYA_SIMPLEX_HPP
#define NOYA_SIMPLEX_HPP 1
/// @complexity Time: Exponential worst case; typically polynomial pivot work per iteration.
/// Space: O(rows * columns).
#include <algorithm>
#include <cassert>
#include <cmath>
#include <limits>
#include <utility>
#include <vector>
namespace noya {
enum class linear_program_status { optimal, infeasible, unbounded };
struct linear_program_result {
linear_program_status sta = linear_program_status::infeasible;
long double obj = 0;
std::vector<long double> sol;
};
namespace simplex_internal {
struct tableau {
static constexpr long double eps = 1e-12L;
// m constraints, n variables.
int m;
int n;
std::vector<int> bas;
std::vector<int> nb;
std::vector<std::vector<long double>> val;
tableau(const std::vector<std::vector<long double>> &a,
const std::vector<long double> &b, const std::vector<long double> &c)
: m(int(b.size())), n(int(c.size())), bas(m), nb(n + 1),
val(m + 2, std::vector<long double>(n + 2)) {
for (int row = 0; row < m; row++) {
for (int col = 0; col < n; col++) {
val[row][col] = a[row][col];
}
bas[row] = n + row;
val[row][n] = -1;
val[row][n + 1] = b[row];
}
for (int col = 0; col < n; col++) {
nb[col] = col;
val[m][col] = -c[col];
}
nb[n] = -1;
val[m + 1][n] = 1;
}
void pivot(int lr, int ec) {
long double inv = 1 / val[lr][ec];
for (int row = 0; row < m + 2; row++) {
if (row == lr) {
continue;
}
for (int col = 0; col < n + 2; col++) {
if (col != ec) {
val[row][col] -= val[lr][col] * val[row][ec] * inv;
}
}
}
for (int col = 0; col < n + 2; col++) {
if (col != ec) {
val[lr][col] *= inv;
}
}
for (int row = 0; row < m + 2; row++) {
if (row != lr) {
val[row][ec] *= -inv;
}
}
val[lr][ec] = inv;
std::swap(bas[lr], nb[ec]);
}
bool optimize(int ph) {
int oi = ph == 1 ? m + 1 : m;
while (true) {
int ent = -1;
for (int col = 0; col <= n; col++) {
if (ph == 2 && nb[col] == -1) {
continue;
}
if (ent == -1 || val[oi][col] < val[oi][ent] - eps ||
(std::abs(val[oi][col] - val[oi][ent]) <= eps &&
nb[col] < nb[ent])) {
ent = col;
}
}
if (ent == -1) {
return true;
}
if (val[oi][ent] >= -eps) {
return true;
}
int lea = -1;
for (int row = 0; row < m; row++) {
if (val[row][ent] <= eps) {
continue;
}
if (lea == -1) {
lea = row;
continue;
}
long double lhs = val[row][n + 1] / val[row][ent];
long double rhs = val[lea][n + 1] / val[lea][ent];
if (lhs < rhs - eps ||
(std::abs(lhs - rhs) <= eps && bas[row] < bas[lea])) {
lea = row;
}
}
if (lea == -1) {
return false;
}
pivot(lea, ent);
}
}
linear_program_result solve() {
int mi = 0;
for (int row = 1; row < m; row++) {
if (val[row][n + 1] < val[mi][n + 1]) {
mi = row;
}
}
if (m > 0 && val[mi][n + 1] < -eps) {
pivot(mi, n);
if (!optimize(1) || val[m + 1][n + 1] < -eps) {
return {linear_program_status::infeasible, 0, {}};
}
if (std::abs(val[m + 1][n + 1]) > eps) {
return {linear_program_status::infeasible, 0, {}};
}
for (int row = 0; row < m; row++) {
if (bas[row] != -1) {
continue;
}
int ent = 0;
for (int col = 1; col <= n; col++) {
if (std::abs(val[row][col]) > std::abs(val[row][ent]) + eps ||
(std::abs(std::abs(val[row][col]) - std::abs(val[row][ent])) <=
eps &&
nb[col] < nb[ent])) {
ent = col;
}
}
if (std::abs(val[row][ent]) > eps) {
pivot(row, ent);
}
}
}
if (!optimize(2)) {
return {linear_program_status::unbounded,
std::numeric_limits<long double>::infinity(),
{}};
}
std::vector<long double> sol(n);
for (int row = 0; row < m; row++) {
if (bas[row] < n) {
sol[bas[row]] = val[row][n + 1];
}
}
return {linear_program_status::optimal, val[m][n + 1], std::move(sol)};
}
};
} // namespace simplex_internal
/// @brief Maximize c*x subject to A*x <= b and x >= 0 with a two-phase
/// simplex tableau; reports optimal, infeasible, or unbounded.
inline linear_program_result
simplex(const std::vector<std::vector<long double>> &a,
const std::vector<long double> &b, const std::vector<long double> &c) {
assert(a.size() == b.size());
for (const auto &row : a) {
assert(row.size() == c.size());
}
return simplex_internal::tableau(a, b, c).solve();
}
} // namespace noya
#endif // NOYA_SIMPLEX_HPP
#include <algorithm>
#include <cassert>
#include <cmath>
#include <limits>
#include <utility>
#include <vector>
/// @complexity Time: Exponential worst case; typically polynomial pivot work per iteration.
/// Space: O(rows * columns).
namespace noya {
enum class linear_program_status { optimal, infeasible, unbounded };
struct linear_program_result {
linear_program_status sta = linear_program_status::infeasible;
long double obj = 0;
std::vector<long double> sol;
};
namespace simplex_internal {
struct tableau {
static constexpr long double eps = 1e-12L;
// m constraints, n variables.
int m;
int n;
std::vector<int> bas;
std::vector<int> nb;
std::vector<std::vector<long double>> val;
tableau(const std::vector<std::vector<long double>> &a,
const std::vector<long double> &b, const std::vector<long double> &c)
: m(int(b.size())), n(int(c.size())), bas(m), nb(n + 1),
val(m + 2, std::vector<long double>(n + 2)) {
for (int row = 0; row < m; row++) {
for (int col = 0; col < n; col++) {
val[row][col] = a[row][col];
}
bas[row] = n + row;
val[row][n] = -1;
val[row][n + 1] = b[row];
}
for (int col = 0; col < n; col++) {
nb[col] = col;
val[m][col] = -c[col];
}
nb[n] = -1;
val[m + 1][n] = 1;
}
void pivot(int lr, int ec) {
long double inv = 1 / val[lr][ec];
for (int row = 0; row < m + 2; row++) {
if (row == lr) {
continue;
}
for (int col = 0; col < n + 2; col++) {
if (col != ec) {
val[row][col] -= val[lr][col] * val[row][ec] * inv;
}
}
}
for (int col = 0; col < n + 2; col++) {
if (col != ec) {
val[lr][col] *= inv;
}
}
for (int row = 0; row < m + 2; row++) {
if (row != lr) {
val[row][ec] *= -inv;
}
}
val[lr][ec] = inv;
std::swap(bas[lr], nb[ec]);
}
bool optimize(int ph) {
int oi = ph == 1 ? m + 1 : m;
while (true) {
int ent = -1;
for (int col = 0; col <= n; col++) {
if (ph == 2 && nb[col] == -1) {
continue;
}
if (ent == -1 || val[oi][col] < val[oi][ent] - eps ||
(std::abs(val[oi][col] - val[oi][ent]) <= eps &&
nb[col] < nb[ent])) {
ent = col;
}
}
if (ent == -1) {
return true;
}
if (val[oi][ent] >= -eps) {
return true;
}
int lea = -1;
for (int row = 0; row < m; row++) {
if (val[row][ent] <= eps) {
continue;
}
if (lea == -1) {
lea = row;
continue;
}
long double lhs = val[row][n + 1] / val[row][ent];
long double rhs = val[lea][n + 1] / val[lea][ent];
if (lhs < rhs - eps ||
(std::abs(lhs - rhs) <= eps && bas[row] < bas[lea])) {
lea = row;
}
}
if (lea == -1) {
return false;
}
pivot(lea, ent);
}
}
linear_program_result solve() {
int mi = 0;
for (int row = 1; row < m; row++) {
if (val[row][n + 1] < val[mi][n + 1]) {
mi = row;
}
}
if (m > 0 && val[mi][n + 1] < -eps) {
pivot(mi, n);
if (!optimize(1) || val[m + 1][n + 1] < -eps) {
return {linear_program_status::infeasible, 0, {}};
}
if (std::abs(val[m + 1][n + 1]) > eps) {
return {linear_program_status::infeasible, 0, {}};
}
for (int row = 0; row < m; row++) {
if (bas[row] != -1) {
continue;
}
int ent = 0;
for (int col = 1; col <= n; col++) {
if (std::abs(val[row][col]) > std::abs(val[row][ent]) + eps ||
(std::abs(std::abs(val[row][col]) - std::abs(val[row][ent])) <=
eps &&
nb[col] < nb[ent])) {
ent = col;
}
}
if (std::abs(val[row][ent]) > eps) {
pivot(row, ent);
}
}
}
if (!optimize(2)) {
return {linear_program_status::unbounded,
std::numeric_limits<long double>::infinity(),
{}};
}
std::vector<long double> sol(n);
for (int row = 0; row < m; row++) {
if (bas[row] < n) {
sol[bas[row]] = val[row][n + 1];
}
}
return {linear_program_status::optimal, val[m][n + 1], std::move(sol)};
}
};
} // namespace simplex_internal
/// @brief Maximize c*x subject to A*x <= b and x >= 0 with a two-phase
/// simplex tableau; reports optimal, infeasible, or unbounded.
inline linear_program_result
simplex(const std::vector<std::vector<long double>> &a,
const std::vector<long double> &b, const std::vector<long double> &c) {
assert(a.size() == b.size());
for (const auto &row : a) {
assert(row.size() == c.size());
}
return simplex_internal::tableau(a, b, c).solve();
}
} // namespace noya