forked from t3nsor/codebook
-
Notifications
You must be signed in to change notification settings - Fork 0
/
simplex.cpp
125 lines (118 loc) · 3.04 KB
/
simplex.cpp
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
// Two-phase simplex algorithm for solving linear programs
// of the form
//
// maximize c^T x
// subject to Ax <= b
// x >= 0
//
// INPUT: A -- an m x n matrix
// b -- an m-dimensional vector
// c -- an n-dimensional vector
// x -- a vector where the optimal solution will be
// stored
//
// OUTPUT: value of the optimal solution (infinity if
// unbounded
// above, nan if infeasible)
//
// To use this code, create an LPSolver object with A, b,
// and c as arguments. Then, call Solve(x).
typedef long double DOUBLE;
typedef vector<DOUBLE> VD;
typedef vector<VD> VVD;
typedef vector<int> VI;
const DOUBLE EPS = 1e-9;
struct LPSolver {
int m, n;
VI B, N;
VVD D;
LPSolver(const VVD &A, const VD &b, const VD &c)
: m(b.size()), n(c.size()), N(n + 1), B(m),
D(m + 2, VD(n + 2)) {
for (int i = 0; i < m; i++)
for (int j = 0; j < n; j++)
D[i][j] = A[i][j];
for (int i = 0; i < m; i++) {
B[i] = n + i;
D[i][n] = -1;
D[i][n + 1] = b[i];
}
for (int j = 0; j < n; j++) {
N[j] = j;
D[m][j] = -c[j];
}
N[n] = -1;
D[m + 1][n] = 1;
}
void Pivot(int r, int s) {
DOUBLE inv = 1.0 / D[r][s];
for (int i = 0; i < m + 2; i++)
if (i != r)
for (int j = 0; j < n + 2; j++)
if (j != s)
D[i][j] -= D[r][j] * D[i][s] * inv;
for (int j = 0; j < n + 2; j++)
if (j != s)
D[r][j] *= inv;
for (int i = 0; i < m + 2; i++)
if (i != r)
D[i][s] *= -inv;
D[r][s] = inv;
swap(B[r], N[s]);
}
bool Simplex(int phase) {
int x = phase == 1 ? m + 1 : m;
while (true) {
int s = -1;
for (int j = 0; j <= n; j++) {
if (phase == 2 && N[j] == -1)
continue;
if (s == -1 || D[x][j] < D[x][s] ||
D[x][j] == D[x][s] && N[j] < N[s])
s = j;
}
if (s < 0 || D[x][s] > -EPS)
return true;
int r = -1;
for (int i = 0; i < m; i++) {
if (D[i][s] < EPS)
continue;
if (r == -1 ||
D[i][n + 1] / D[i][s] < D[r][n + 1] / D[r][s] ||
D[i][n + 1] / D[i][s] == D[r][n + 1] / D[r][s] &&
B[i] < B[r])
r = i;
}
if (r == -1)
return false;
Pivot(r, s);
}
}
DOUBLE Solve(VD &x) {
int r = 0;
for (int i = 1; i < m; i++)
if (D[i][n + 1] < D[r][n + 1])
r = i;
if (D[r][n + 1] <= -EPS) {
Pivot(r, n);
if (!Simplex(1) || D[m + 1][n + 1] < -EPS)
return -numeric_limits<DOUBLE>::infinity();
for (int i = 0; i < m; i++)
if (B[i] == -1) {
int s = -1;
for (int j = 0; j <= n; j++)
if (s == -1 || D[i][j] < D[i][s] ||
D[i][j] == D[i][s] && N[j] < N[s])
s = j;
Pivot(i, s);
}
}
if (!Simplex(2))
return numeric_limits<DOUBLE>::infinity();
x = VD(n);
for (int i = 0; i < m; i++)
if (B[i] < n)
x[B[i]] = D[i][n + 1];
return D[m][n + 1];
}
};