-
Notifications
You must be signed in to change notification settings - Fork 1
/
test.py
60 lines (52 loc) · 2.05 KB
/
test.py
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
import numpy as np
import cvxpy as cp
from lp_solver import LPSolver
def generate_data(m, n, feasible=True):
"""
Randomly generates an LP problem
Returns A (m x n), b (m), c (n)
"""
assert m < n
A = np.random.randn(m, n)
A[0,:] = np.random.rand(n) + 0.1 # make sure problem is bounded
if feasible:
p = np.random.rand(n) + 0.01 # all positive
else:
p = np.random.randn(n)
b = np.dot(A, p)
c = np.random.rand(n)
return A, b, c
if __name__ == "__main__":
# Generate a feasible LP problem
m, n = 100, 500
A, b, c = generate_data(m, n, feasible=True)
# Solve the problem using LPSolver
solver = LPSolver(mu=10, tol=1e-4)
solver.solve(A, b, c)
assert solver.status == 'optimal'
x_opt = solver.x_opt
print("LP Solver optimal value: {}".format(solver.value))
print("LPSolver number of centering steps: {}".format(solver.num_steps))
print("LPSolver norm(Ax-b): {}".format(np.linalg.norm(A.dot(x_opt) - b)))
print("LPSolver number of negative x_i: {}".format((x_opt < 0).sum()))
# Solve the same problem using CVXPY, and compare answers
x_cp = cp.Variable(n)
objective = c.T * x_cp
problem = cp.Problem(cp.Minimize(objective), [A * x_cp == b, x_cp >= 0])
problem.solve()
assert problem.status == 'optimal'
print("\nCVX optimal value: {}".format(objective.value))
print("Percent diff: {}%".format((solver.value - objective.value) / objective.value * 100))
# Experiment infeasible problems
num_runs = 10
print("\nTest LPSolver on {} infeasible problems...".format(num_runs))
for i in range(num_runs):
A, b, c = generate_data(m, n, feasible=False)
solver = LPSolver(mu=10, tol=1e-4)
solver.solve(A, b, c)
problem = cp.Problem(cp.Minimize(objective), [A * x_cp == b, x_cp >= 0])
problem.solve()
if problem.status != 'infeasible':
print("WHAT?!")
assert solver.status == 'infeasible'
print("LPSolver and CVXPY are consistent on the {} infeasible problems.".format(num_runs))