Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

TST: optimize.linprog/milp: add tests for various bug reports #20589

Draft
wants to merge 1 commit into
base: main
Choose a base branch
from
Draft
Show file tree
Hide file tree
Changes from all commits
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Jump to
Jump to file
Failed to load files.
Diff view
Diff view
95 changes: 81 additions & 14 deletions scipy/optimize/tests/test_linprog.py
Original file line number Diff line number Diff line change
Expand Up @@ -1703,20 +1703,6 @@ def test_bug_10466(self):
method=self.method, options=o)
assert_allclose(res.fun, -8589934560)

def test_bug_20584(self):
"""
Test that when integrality is a list of all zeros, linprog gives the
same result as when it is an array of all zeros / integrality=None
"""
c = [1, 1]
A_ub = [[-1, 0]]
b_ub = [-2.5]
res1 = linprog(c, A_ub=A_ub, b_ub=b_ub, integrality=[0, 0])
res2 = linprog(c, A_ub=A_ub, b_ub=b_ub, integrality=np.asarray([0, 0]))
res3 = linprog(c, A_ub=A_ub, b_ub=b_ub, integrality=None)
assert_equal(res1.x, res2.x)
assert_equal(res1.x, res3.x)


#########################
# Method-specific Tests #
Expand Down Expand Up @@ -1877,6 +1863,73 @@ def test_complementary_slackness(self):
# non-zero RHS
assert np.allclose(res.ineqlin.marginals @ (b_ub - A_ub @ res.x), 0)

def test_bug_20336(self):
"""
Test that `linprog` now solves a poorly-scaled problem
"""
boundaries = [(10000.0, 9010000.0), (0.0, None), (10000.0, None),
(0.0, 84.62623413258109), (10000.0, None), (10000.0, None),
(10000.0, None), (10000.0, None), (10000.0, None),
(10000.0, None), (10000.0, None), (10000.0, None),
(10000.0, None), (None, None), (None, None), (None, None),
(None, None), (None, None), (None, None), (None, None),
(None, None), (None, None), (None, None), (None, None),
(None, None), (None, None), (None, None), (None, None),
(None, None), (None, None), (None, None), (None, None),
(None, None)]
eq_entries = [-1.0, 1.0, -1.0, 1.0, -1.0, 1.0, -1.0, 1.0, -1.0, 1.0, -1.0, 1.0,
-1.0, 1.0, -1.0, 1.0, -1.0, 1.0, -1.0, 1.0, 1.0, 1.0, -1.0, 0.001,
-0.001, 3.7337777768059636e-10, 3.7337777768059636e-10, 1.0, -1.0,
0.001, -0.001, 3.7337777768059636e-10, 3.7337777768059636e-10,
1.0, -1.0, 0.001, -0.001, 3.7337777768059636e-10,
3.7337777768059636e-10, 1.0, -1.0, 0.001, -0.001,
3.7337777768059636e-10, 3.7337777768059636e-10, 1.0, -1.0, 0.001,
-0.001, 3.7337777768059636e-10, 3.7337777768059636e-10, 1.0, -1.0,
0.001, -0.001, 3.7337777768059636e-10, 3.7337777768059636e-10,
1.0, -1.0, 0.001, -0.001, 3.7337777768059636e-10,
3.7337777768059636e-10, 1.0, -1.0, 0.001, -0.001,
3.7337777768059636e-10, 3.7337777768059636e-10, 1.0, -1.0, 0.001,
-0.001, 3.7337777768059636e-10, 3.7337777768059636e-10, 1.0,
-1.0, 0.001, -0.001, 3.7337777768059636e-10,
3.7337777768059636e-10, 1.0, -1.0]
eq_indizes = [0, 0, 1, 1, 2, 2, 3, 3, 4, 4, 5, 5, 6, 6, 7, 7, 8, 8, 9, 9, 10,
11, 11, 12, 12, 12, 12, 13, 13, 14, 14, 14, 14, 15, 15, 16, 16,
16, 16, 17, 17, 18, 18, 18, 18, 19, 19, 20, 20, 20, 20, 21, 21,
22, 22, 22, 22, 23, 23, 24, 24, 24, 24, 25, 25, 26, 26, 26, 26,
27, 27, 28, 28, 28, 28, 29, 29, 30, 30, 30, 30, 31, 31]
eq_vars = [15, 14, 17, 16, 19, 18, 21, 20, 23, 22, 25, 24, 27, 26, 29, 28, 31,
30, 13, 1, 0, 32, 3, 14, 13, 4, 0, 4, 0, 32, 31, 2, 12, 2, 12, 16,
15, 5, 4, 5, 4, 18, 17, 6, 5, 6, 5, 20, 19, 7, 6, 7, 6, 22, 21, 8,
7, 8, 7, 24, 23, 9, 8, 9, 8, 26, 25, 10, 9, 10, 9, 28, 27, 11, 10,
11, 10, 30, 29, 12, 11, 12, 11]
eq_values = [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 9000000.0, 0.0,
0.006587392118285457, -5032.197406716549, 0.0041860502789104696,
-7918.93439542944, 0.0063205763583549035, -5244.625751707402,
0.006053760598424349, -5475.7793929428, 0.005786944838493795,
-5728.248403917573, 0.0055201290785632405, -6005.123623538355,
0.005253313318632687, -6310.123825488683, 0.004986497558702133,
-6647.763714796453, 0.004719681798771578, -7023.578908071522,
0.004452866038841024, -7444.431798646482]
coefficients = [0.0, 0.0, 0.0, -0.011816666666666668, 0.0, 0.0, 0.0, 0.0, 0.0,
0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0,
0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0]
np_eq_entries = np.asarray(eq_entries, dtype=np.float64)
np_eq_indizes = np.asarray(eq_indizes, dtype=np.int32)
np_eq_vars = np.asarray(eq_vars, dtype=np.int32)

a_eq= scipy.sparse.csr_array((np_eq_entries,(np_eq_indizes, np_eq_vars)),
shape=(32, 33))
b_eq = np.asarray(eq_values, dtype=np.float64)
c = np.asarray(coefficients, dtype=np.float64)

result = scipy.optimize.linprog(c, A_ub=None, b_ub=None, A_eq=a_eq, b_eq=b_eq,
bounds=boundaries)
assert result.status==0
x = result.x
n_r_x = np.linalg.norm(a_eq @ x - b_eq)
n_r = np.linalg.norm(result.eqlin.residual)
assert_allclose(n_r, n_r_x)


################################
# Simplex Option-Specific Tests#
Expand Down Expand Up @@ -2423,6 +2476,20 @@ def test_semi_continuous(self):
np.testing.assert_allclose(res.x, [0, 0, 1.5, 1])
assert res.status == 0

def test_bug_20584(self):
"""
Test that when integrality is a list of all zeros, linprog gives the
same result as when it is an array of all zeros / integrality=None
"""
c = [1, 1]
A_ub = [[-1, 0]]
b_ub = [-2.5]
res1 = linprog(c, A_ub=A_ub, b_ub=b_ub, integrality=[0, 0])
res2 = linprog(c, A_ub=A_ub, b_ub=b_ub, integrality=np.asarray([0, 0]))
res3 = linprog(c, A_ub=A_ub, b_ub=b_ub, integrality=None)
assert_equal(res1.x, res2.x)
assert_equal(res1.x, res3.x)


###########################
# Autoscale-Specific Tests#
Expand Down
69 changes: 69 additions & 0 deletions scipy/optimize/tests/test_milp.py
Original file line number Diff line number Diff line change
Expand Up @@ -385,3 +385,72 @@ def test_mip_rel_gap_passdown():
# check that differences between solution gaps are declining
# monotonically with the mip_rel_gap parameter.
assert np.all(np.diff(sol_mip_gaps) < 0)

def test_large_numbers_gh20116():
h = 10 ** 12
A = np.array([[100.4534, h], [100.4534, -h]])
b = np.array([h, 0])
constraints = LinearConstraint(A=A, ub=b)
bounds = Bounds([0, 0], [1, 1])
c = np.array([0, 0])
res = milp(c=c, constraints=constraints, bounds=bounds, integrality=1)
assert res.status == 0
assert np.all(A @ res.x < b)


def test_presolve_gh18907():
from scipy.optimize import milp
import numpy as np
inf = np.inf

# set up problem
c = np.array([-0.85850509, -0.82892676, -0.80026454, -0.63015535, -0.5099006,
-0.50077193, -0.4894404, -0.47285865, -0.39867774, -0.38069646,
-0.36733012, -0.36733012, -0.35820411, -0.31576141, -0.20626091,
-0.12466144, -0.10679516, -0.1061887, -0.1061887, -0.1061887,
-0., -0., -0., -0., 0., 0., 0., 0.])

A = np.array([[1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1.,
1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 0., 0., 0., 0.],
[0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.,
1., 0., 0., 0., 0., 0., 1., 0., 0., 0., -25., -0., -0., -0.],
[0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.,
-1., 0., 0., 0., 0., 0., -1., 0., 0., 0., 2., 0., 0., 0.],
[0., 0., 0., 0., 1., 1., 1., 1., 0., 1., 0., 0., 0., 0., 0.,
0., 0., 0., 0., 0., 0., 0., 0., 0., -0., -25., -0., -0.],
[0., 0., 0., 0., -1., -1., -1., -1., 0., -1., 0., 0., 0.,
0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 2., 0., 0.],
[0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.,
0., 0., 1., 1., 1., 0., 0., 0., 0., -0., -0., -25., -0.],
[0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.,
0., 0., -1., -1., -1., 0., 0., 0., 0., 0., 0., 2., 0.],
[1., 1., 1., 1., 0., 0., 0., 0., 1., 0., 1., 1., 1., 1., 0.,
1., 1., 0., 0., 0., 0., 1., 1., 1., -0., -0., -0., -25.],
[-1., -1., -1., -1., 0., 0., 0., 0., -1., 0., -1., -1., -1., -1.,
0., -1., -1., 0., 0., 0., 0., -1., -1., -1., 0., 0., 0., 2.]])
bl = np.array([-inf, -inf, -inf, -inf, -inf, -inf, -inf, -inf, -inf])
bu = np.array([100., 0., 0., 0., 0., 0., 0., 0., 0.])
constraints = LinearConstraint(A, bl, bu)
integrality = 1
bounds = (0, 1)
r1 = milp(c=c, constraints=constraints, integrality=integrality, bounds=bounds,
options={'presolve': True})
r2 = milp(c=c, constraints=constraints, integrality=integrality, bounds=bounds,
options={'presolve': False})
assert r1.status == r2.status
assert_allclose(r1.x, r2.x)

# another example from the same issue
bounds = Bounds(lb=0, ub=1)
integrality = [1, 1, 0, 0]
c = [10, 9.52380952, -1000, -952.38095238]
A = [[1, 1, 0, 0], [0, 0, 1, 1], [200, 0, 0, 0], [0, 200, 0, 0],
[0, 0, 2000, 0], [0, 0, 0, 2000], [-1, 0, 1, 0], [-1, -1, 0, 1]]
ub = [1, 1, 200, 200, 1000, 1000, 0, 0]
constraints = LinearConstraint(A, ub=ub)
r1 = milp(c=c, constraints=constraints, bounds=bounds,
integrality=integrality, options={"presolve": False})
r2 = milp(c=c, constraints=constraints, bounds=bounds,
integrality=integrality, options={"presolve": False})
assert r1.status == r2.status
assert_allclose(r1.x, r2.x)