test_linprog.py 53 KB

1234567891011121314151617181920212223242526272829303132333435363738394041424344454647484950515253545556575859606162636465666768697071727374757677787980818283848586878889909192939495969798991001011021031041051061071081091101111121131141151161171181191201211221231241251261271281291301311321331341351361371381391401411421431441451461471481491501511521531541551561571581591601611621631641651661671681691701711721731741751761771781791801811821831841851861871881891901911921931941951961971981992002012022032042052062072082092102112122132142152162172182192202212222232242252262272282292302312322332342352362372382392402412422432442452462472482492502512522532542552562572582592602612622632642652662672682692702712722732742752762772782792802812822832842852862872882892902912922932942952962972982993003013023033043053063073083093103113123133143153163173183193203213223233243253263273283293303313323333343353363373383393403413423433443453463473483493503513523533543553563573583593603613623633643653663673683693703713723733743753763773783793803813823833843853863873883893903913923933943953963973983994004014024034044054064074084094104114124134144154164174184194204214224234244254264274284294304314324334344354364374384394404414424434444454464474484494504514524534544554564574584594604614624634644654664674684694704714724734744754764774784794804814824834844854864874884894904914924934944954964974984995005015025035045055065075085095105115125135145155165175185195205215225235245255265275285295305315325335345355365375385395405415425435445455465475485495505515525535545555565575585595605615625635645655665675685695705715725735745755765775785795805815825835845855865875885895905915925935945955965975985996006016026036046056066076086096106116126136146156166176186196206216226236246256266276286296306316326336346356366376386396406416426436446456466476486496506516526536546556566576586596606616626636646656666676686696706716726736746756766776786796806816826836846856866876886896906916926936946956966976986997007017027037047057067077087097107117127137147157167177187197207217227237247257267277287297307317327337347357367377387397407417427437447457467477487497507517527537547557567577587597607617627637647657667677687697707717727737747757767777787797807817827837847857867877887897907917927937947957967977987998008018028038048058068078088098108118128138148158168178188198208218228238248258268278288298308318328338348358368378388398408418428438448458468478488498508518528538548558568578588598608618628638648658668678688698708718728738748758768778788798808818828838848858868878888898908918928938948958968978988999009019029039049059069079089099109119129139149159169179189199209219229239249259269279289299309319329339349359369379389399409419429439449459469479489499509519529539549559569579589599609619629639649659669679689699709719729739749759769779789799809819829839849859869879889899909919929939949959969979989991000100110021003100410051006100710081009101010111012101310141015101610171018101910201021102210231024102510261027102810291030103110321033103410351036103710381039104010411042104310441045104610471048104910501051105210531054105510561057105810591060106110621063106410651066106710681069107010711072107310741075107610771078107910801081108210831084108510861087108810891090109110921093109410951096109710981099110011011102110311041105110611071108110911101111111211131114111511161117111811191120112111221123112411251126112711281129113011311132113311341135113611371138113911401141114211431144114511461147114811491150115111521153115411551156115711581159116011611162116311641165116611671168116911701171117211731174117511761177117811791180118111821183118411851186118711881189119011911192119311941195119611971198119912001201120212031204120512061207120812091210121112121213121412151216121712181219122012211222122312241225122612271228122912301231123212331234123512361237123812391240124112421243124412451246124712481249125012511252125312541255125612571258125912601261126212631264126512661267126812691270127112721273127412751276127712781279128012811282128312841285128612871288128912901291129212931294129512961297129812991300130113021303130413051306130713081309131013111312131313141315131613171318131913201321132213231324132513261327132813291330133113321333133413351336133713381339134013411342134313441345134613471348134913501351135213531354135513561357135813591360136113621363136413651366136713681369137013711372137313741375
  1. """
  2. Unit test for Linear Programming
  3. """
  4. from __future__ import division, print_function, absolute_import
  5. import numpy as np
  6. from numpy.testing import (assert_, assert_allclose, assert_equal,
  7. assert_array_less)
  8. from pytest import raises as assert_raises
  9. from scipy.optimize import linprog, OptimizeWarning
  10. from scipy._lib._numpy_compat import _assert_warns, suppress_warnings
  11. from scipy.sparse.linalg import MatrixRankWarning
  12. import pytest
  13. def magic_square(n):
  14. np.random.seed(0)
  15. M = n * (n**2 + 1) / 2
  16. numbers = np.arange(n**4) // n**2 + 1
  17. numbers = numbers.reshape(n**2, n, n)
  18. zeros = np.zeros((n**2, n, n))
  19. A_list = []
  20. b_list = []
  21. # Rule 1: use every number exactly once
  22. for i in range(n**2):
  23. A_row = zeros.copy()
  24. A_row[i, :, :] = 1
  25. A_list.append(A_row.flatten())
  26. b_list.append(1)
  27. # Rule 2: Only one number per square
  28. for i in range(n):
  29. for j in range(n):
  30. A_row = zeros.copy()
  31. A_row[:, i, j] = 1
  32. A_list.append(A_row.flatten())
  33. b_list.append(1)
  34. # Rule 3: sum of rows is M
  35. for i in range(n):
  36. A_row = zeros.copy()
  37. A_row[:, i, :] = numbers[:, i, :]
  38. A_list.append(A_row.flatten())
  39. b_list.append(M)
  40. # Rule 4: sum of columns is M
  41. for i in range(n):
  42. A_row = zeros.copy()
  43. A_row[:, :, i] = numbers[:, :, i]
  44. A_list.append(A_row.flatten())
  45. b_list.append(M)
  46. # Rule 5: sum of diagonals is M
  47. A_row = zeros.copy()
  48. A_row[:, range(n), range(n)] = numbers[:, range(n), range(n)]
  49. A_list.append(A_row.flatten())
  50. b_list.append(M)
  51. A_row = zeros.copy()
  52. A_row[:, range(n), range(-1, -n - 1, -1)] = \
  53. numbers[:, range(n), range(-1, -n - 1, -1)]
  54. A_list.append(A_row.flatten())
  55. b_list.append(M)
  56. A = np.array(np.vstack(A_list), dtype=float)
  57. b = np.array(b_list, dtype=float)
  58. c = np.random.rand(A.shape[1])
  59. return A, b, c, numbers
  60. def lpgen_2d(m, n):
  61. """ -> A b c LP test: m*n vars, m+n constraints
  62. row sums == n/m, col sums == 1
  63. https://gist.github.com/denis-bz/8647461
  64. """
  65. np.random.seed(0)
  66. c = - np.random.exponential(size=(m, n))
  67. Arow = np.zeros((m, m * n))
  68. brow = np.zeros(m)
  69. for j in range(m):
  70. j1 = j + 1
  71. Arow[j, j * n:j1 * n] = 1
  72. brow[j] = n / m
  73. Acol = np.zeros((n, m * n))
  74. bcol = np.zeros(n)
  75. for j in range(n):
  76. j1 = j + 1
  77. Acol[j, j::n] = 1
  78. bcol[j] = 1
  79. A = np.vstack((Arow, Acol))
  80. b = np.hstack((brow, bcol))
  81. return A, b, c.ravel()
  82. def _assert_iteration_limit_reached(res, maxiter):
  83. assert_(not res.success, "Incorrectly reported success")
  84. assert_(res.success < maxiter, "Incorrectly reported number of iterations")
  85. assert_equal(res.status, 1, "Failed to report iteration limit reached")
  86. def _assert_infeasible(res):
  87. # res: linprog result object
  88. assert_(not res.success, "incorrectly reported success")
  89. assert_equal(res.status, 2, "failed to report infeasible status")
  90. def _assert_unbounded(res):
  91. # res: linprog result object
  92. assert_(not res.success, "incorrectly reported success")
  93. assert_equal(res.status, 3, "failed to report unbounded status")
  94. def _assert_unable_to_find_basic_feasible_sol(res):
  95. # res: linprog result object
  96. assert_(not res.success, "incorrectly reported success")
  97. assert_equal(res.status, 2, "failed to report optimization failure")
  98. def _assert_success(res, desired_fun=None, desired_x=None,
  99. rtol=1e-8, atol=1e-8):
  100. # res: linprog result object
  101. # desired_fun: desired objective function value or None
  102. # desired_x: desired solution or None
  103. if not res.success:
  104. msg = "linprog status {0}, message: {1}".format(res.status,
  105. res.message)
  106. raise AssertionError(msg)
  107. assert_equal(res.status, 0)
  108. if desired_fun is not None:
  109. assert_allclose(res.fun, desired_fun,
  110. err_msg="converged to an unexpected objective value",
  111. rtol=rtol, atol=atol)
  112. if desired_x is not None:
  113. assert_allclose(res.x, desired_x,
  114. err_msg="converged to an unexpected solution",
  115. rtol=rtol, atol=atol)
  116. class LinprogCommonTests(object):
  117. def test_docstring_example(self):
  118. # Example from linprog docstring.
  119. c = [-1, 4]
  120. A = [[-3, 1], [1, 2]]
  121. b = [6, 4]
  122. x0_bounds = (None, None)
  123. x1_bounds = (-3, None)
  124. res = linprog(c, A_ub=A, b_ub=b, bounds=(x0_bounds, x1_bounds),
  125. options=self.options, method=self.method)
  126. _assert_success(res, desired_fun=-22)
  127. def test_aliasing_b_ub(self):
  128. c = np.array([1.0])
  129. A_ub = np.array([[1.0]])
  130. b_ub_orig = np.array([3.0])
  131. b_ub = b_ub_orig.copy()
  132. bounds = (-4.0, np.inf)
  133. res = linprog(c, A_ub=A_ub, b_ub=b_ub, bounds=bounds,
  134. method=self.method, options=self.options)
  135. _assert_success(res, desired_fun=-4, desired_x=[-4])
  136. assert_allclose(b_ub_orig, b_ub)
  137. def test_aliasing_b_eq(self):
  138. c = np.array([1.0])
  139. A_eq = np.array([[1.0]])
  140. b_eq_orig = np.array([3.0])
  141. b_eq = b_eq_orig.copy()
  142. bounds = (-4.0, np.inf)
  143. res = linprog(c, A_eq=A_eq, b_eq=b_eq, bounds=bounds,
  144. method=self.method, options=self.options)
  145. _assert_success(res, desired_fun=3, desired_x=[3])
  146. assert_allclose(b_eq_orig, b_eq)
  147. def test_bounds_second_form_unbounded_below(self):
  148. c = np.array([1.0])
  149. A_eq = np.array([[1.0]])
  150. b_eq = np.array([3.0])
  151. bounds = (None, 10.0)
  152. res = linprog(c, A_eq=A_eq, b_eq=b_eq, bounds=bounds,
  153. method=self.method, options=self.options)
  154. _assert_success(res, desired_fun=3, desired_x=[3])
  155. def test_bounds_second_form_unbounded_above(self):
  156. c = np.array([1.0])
  157. A_eq = np.array([[1.0]])
  158. b_eq = np.array([3.0])
  159. bounds = (1.0, None)
  160. res = linprog(c, A_eq=A_eq, b_eq=b_eq, bounds=bounds,
  161. method=self.method, options=self.options)
  162. _assert_success(res, desired_fun=3, desired_x=[3])
  163. def test_non_ndarray_args(self):
  164. c = [1.0]
  165. A_ub = [[1.0]]
  166. b_ub = [3.0]
  167. A_eq = [[1.0]]
  168. b_eq = [2.0]
  169. bounds = (-1.0, 10.0)
  170. res = linprog(c, A_ub=A_ub, b_ub=b_ub, A_eq=A_eq, b_eq=b_eq,
  171. bounds=bounds, method=self.method, options=self.options)
  172. _assert_success(res, desired_fun=2, desired_x=[2])
  173. def test_linprog_upper_bound_constraints(self):
  174. # Maximize a linear function subject to only linear upper bound
  175. # constraints.
  176. # http://www.dam.brown.edu/people/huiwang/classes/am121/Archive/simplex_121_c.pdf
  177. c = np.array([3, 2]) * -1 # maximize
  178. A_ub = [[2, 1],
  179. [1, 1],
  180. [1, 0]]
  181. b_ub = [10, 8, 4]
  182. res = (linprog(c, A_ub=A_ub, b_ub=b_ub,
  183. method=self.method, options=self.options))
  184. _assert_success(res, desired_fun=-18, desired_x=[2, 6])
  185. def test_linprog_mixed_constraints(self):
  186. # Minimize linear function subject to non-negative variables.
  187. # http://www.statslab.cam.ac.uk/~ff271/teaching/opt/notes/notes8.pdf
  188. # (dead link)
  189. c = [6, 3]
  190. A_ub = [[0, 3],
  191. [-1, -1],
  192. [-2, 1]]
  193. b_ub = [2, -1, -1]
  194. res = linprog(c, A_ub=A_ub, b_ub=b_ub,
  195. method=self.method, options=self.options)
  196. _assert_success(res, desired_fun=5, desired_x=[2 / 3, 1 / 3])
  197. def test_linprog_cyclic_recovery(self):
  198. # Test linprogs recovery from cycling using the Klee-Minty problem
  199. # Klee-Minty https://www.math.ubc.ca/~israel/m340/kleemin3.pdf
  200. c = np.array([100, 10, 1]) * -1 # maximize
  201. A_ub = [[1, 0, 0],
  202. [20, 1, 0],
  203. [200, 20, 1]]
  204. b_ub = [1, 100, 10000]
  205. res = linprog(c, A_ub=A_ub, b_ub=b_ub,
  206. method=self.method, options=self.options)
  207. _assert_success(res, desired_x=[0, 0, 10000], atol=5e-6, rtol=1e-7)
  208. def test_linprog_cyclic_bland(self):
  209. # Test the effect of Bland's rule on a cycling problem
  210. c = np.array([-10, 57, 9, 24.])
  211. A_ub = np.array([[0.5, -5.5, -2.5, 9],
  212. [0.5, -1.5, -0.5, 1],
  213. [1, 0, 0, 0]])
  214. b_ub = [0, 0, 1]
  215. maxiter = 100
  216. o = {key: val for key, val in self.options.items()}
  217. o['maxiter'] = maxiter
  218. res = linprog(c, A_ub=A_ub, b_ub=b_ub, options=o,
  219. method=self.method)
  220. if self.method == 'simplex'and not self.options.get('bland'):
  221. _assert_iteration_limit_reached(res, o['maxiter'])
  222. else:
  223. _assert_success(res, desired_x=[1, 0, 1, 0])
  224. def test_linprog_cyclic_bland_bug_8561(self):
  225. # Test that pivot row is chosen correctly when using Bland's rule
  226. c = np.array([7, 0, -4, 1.5, 1.5])
  227. A_ub = np.array([
  228. [4, 5.5, 1.5, 1.0, -3.5],
  229. [1, -2.5, -2, 2.5, 0.5],
  230. [3, -0.5, 4, -12.5, -7],
  231. [-1, 4.5, 2, -3.5, -2],
  232. [5.5, 2, -4.5, -1, 9.5]])
  233. b_ub = np.array([0, 0, 0, 0, 1])
  234. if self.method == "simplex":
  235. res = linprog(c, A_ub=A_ub, b_ub=b_ub,
  236. options=dict(maxiter=100, bland=True),
  237. method=self.method)
  238. else:
  239. res = linprog(c, A_ub=A_ub, b_ub=b_ub, options=self.options,
  240. method=self.method)
  241. _assert_success(res, desired_x=[0, 0, 19, 16/3, 29/3])
  242. def test_linprog_unbounded(self):
  243. # Test linprog response to an unbounded problem
  244. c = np.array([1, 1]) * -1 # maximize
  245. A_ub = [[-1, 1],
  246. [-1, -1]]
  247. b_ub = [-1, -2]
  248. res = linprog(c, A_ub=A_ub, b_ub=b_ub,
  249. method=self.method, options=self.options)
  250. _assert_unbounded(res)
  251. def test_linprog_infeasible(self):
  252. # Test linrpog response to an infeasible problem
  253. c = [-1, -1]
  254. A_ub = [[1, 0],
  255. [0, 1],
  256. [-1, -1]]
  257. b_ub = [2, 2, -5]
  258. res = linprog(c, A_ub=A_ub, b_ub=b_ub,
  259. method=self.method, options=self.options)
  260. _assert_infeasible(res)
  261. def test_nontrivial_problem(self):
  262. # Test linprog for a problem involving all constraint types,
  263. # negative resource limits, and rounding issues.
  264. c = [-1, 8, 4, -6]
  265. A_ub = [[-7, -7, 6, 9],
  266. [1, -1, -3, 0],
  267. [10, -10, -7, 7],
  268. [6, -1, 3, 4]]
  269. b_ub = [-3, 6, -6, 6]
  270. A_eq = [[-10, 1, 1, -8]]
  271. b_eq = [-4]
  272. res = linprog(c, A_ub=A_ub, b_ub=b_ub, A_eq=A_eq, b_eq=b_eq,
  273. method=self.method, options=self.options)
  274. _assert_success(res, desired_fun=7083 / 1391,
  275. desired_x=[101 / 1391, 1462 / 1391, 0, 752 / 1391])
  276. def test_negative_variable(self):
  277. # Test linprog with a problem with one unbounded variable and
  278. # another with a negative lower bound.
  279. c = np.array([-1, 4]) * -1 # maximize
  280. A_ub = np.array([[-3, 1],
  281. [1, 2]], dtype=np.float64)
  282. A_ub_orig = A_ub.copy()
  283. b_ub = [6, 4]
  284. x0_bounds = (-np.inf, np.inf)
  285. x1_bounds = (-3, np.inf)
  286. res = linprog(c, A_ub=A_ub, b_ub=b_ub, bounds=(x0_bounds, x1_bounds),
  287. method=self.method, options=self.options)
  288. assert_equal(A_ub, A_ub_orig) # user input not overwritten
  289. _assert_success(res, desired_fun=-80 / 7, desired_x=[-8 / 7, 18 / 7])
  290. def test_large_problem(self):
  291. # Test linprog simplex with a rather large problem (400 variables,
  292. # 40 constraints) generated by https://gist.github.com/denis-bz/8647461
  293. A, b, c = lpgen_2d(20, 20)
  294. res = linprog(c, A_ub=A, b_ub=b,
  295. method=self.method, options=self.options)
  296. _assert_success(res, desired_fun=-64.049494229)
  297. def test_network_flow(self):
  298. # A network flow problem with supply and demand at nodes
  299. # and with costs along directed edges.
  300. # https://www.princeton.edu/~rvdb/542/lectures/lec10.pdf
  301. c = [2, 4, 9, 11, 4, 3, 8, 7, 0, 15, 16, 18]
  302. n, p = -1, 1
  303. A_eq = [
  304. [n, n, p, 0, p, 0, 0, 0, 0, p, 0, 0],
  305. [p, 0, 0, p, 0, p, 0, 0, 0, 0, 0, 0],
  306. [0, 0, n, n, 0, 0, 0, 0, 0, 0, 0, 0],
  307. [0, 0, 0, 0, 0, 0, p, p, 0, 0, p, 0],
  308. [0, 0, 0, 0, n, n, n, 0, p, 0, 0, 0],
  309. [0, 0, 0, 0, 0, 0, 0, n, n, 0, 0, p],
  310. [0, 0, 0, 0, 0, 0, 0, 0, 0, n, n, n]]
  311. b_eq = [0, 19, -16, 33, 0, 0, -36]
  312. res = linprog(c=c, A_eq=A_eq, b_eq=b_eq,
  313. method=self.method, options=self.options)
  314. _assert_success(res, desired_fun=755, atol=1e-6, rtol=1e-7)
  315. def test_network_flow_limited_capacity(self):
  316. # A network flow problem with supply and demand at nodes
  317. # and with costs and capacities along directed edges.
  318. # http://blog.sommer-forst.de/2013/04/10/
  319. cost = [2, 2, 1, 3, 1]
  320. bounds = [
  321. [0, 4],
  322. [0, 2],
  323. [0, 2],
  324. [0, 3],
  325. [0, 5]]
  326. n, p = -1, 1
  327. A_eq = [
  328. [n, n, 0, 0, 0],
  329. [p, 0, n, n, 0],
  330. [0, p, p, 0, n],
  331. [0, 0, 0, p, p]]
  332. b_eq = [-4, 0, 0, 4]
  333. with suppress_warnings() as sup:
  334. sup.filter(RuntimeWarning, "scipy.linalg.solve\nIll...")
  335. sup.filter(OptimizeWarning, "A_eq does not appear...")
  336. sup.filter(OptimizeWarning, "Solving system with option...")
  337. res = linprog(c=cost, A_eq=A_eq, b_eq=b_eq, bounds=bounds,
  338. method=self.method, options=self.options)
  339. _assert_success(res, desired_fun=14)
  340. def test_simplex_algorithm_wikipedia_example(self):
  341. # https://en.wikipedia.org/wiki/Simplex_algorithm#Example
  342. Z = [-2, -3, -4]
  343. A_ub = [
  344. [3, 2, 1],
  345. [2, 5, 3]]
  346. b_ub = [10, 15]
  347. res = linprog(c=Z, A_ub=A_ub, b_ub=b_ub,
  348. method=self.method, options=self.options)
  349. _assert_success(res, desired_fun=-20)
  350. def test_enzo_example(self):
  351. # https://github.com/scipy/scipy/issues/1779 lp2.py
  352. #
  353. # Translated from Octave code at:
  354. # http://www.ecs.shimane-u.ac.jp/~kyoshida/lpeng.htm
  355. # and placed under MIT licence by Enzo Michelangeli
  356. # with permission explicitly granted by the original author,
  357. # Prof. Kazunobu Yoshida
  358. c = [4, 8, 3, 0, 0, 0]
  359. A_eq = [
  360. [2, 5, 3, -1, 0, 0],
  361. [3, 2.5, 8, 0, -1, 0],
  362. [8, 10, 4, 0, 0, -1]]
  363. b_eq = [185, 155, 600]
  364. res = linprog(c=c, A_eq=A_eq, b_eq=b_eq,
  365. method=self.method, options=self.options)
  366. _assert_success(res, desired_fun=317.5,
  367. desired_x=[66.25, 0, 17.5, 0, 183.75, 0],
  368. atol=6e-6, rtol=1e-7)
  369. def test_enzo_example_b(self):
  370. # rescued from https://github.com/scipy/scipy/pull/218
  371. c = [2.8, 6.3, 10.8, -2.8, -6.3, -10.8]
  372. A_eq = [[-1, -1, -1, 0, 0, 0],
  373. [0, 0, 0, 1, 1, 1],
  374. [1, 0, 0, 1, 0, 0],
  375. [0, 1, 0, 0, 1, 0],
  376. [0, 0, 1, 0, 0, 1]]
  377. b_eq = [-0.5, 0.4, 0.3, 0.3, 0.3]
  378. with suppress_warnings() as sup:
  379. sup.filter(OptimizeWarning, "A_eq does not appear...")
  380. res = linprog(c=c, A_eq=A_eq, b_eq=b_eq,
  381. method=self.method, options=self.options)
  382. _assert_success(res, desired_fun=-1.77,
  383. desired_x=[0.3, 0.2, 0.0, 0.0, 0.1, 0.3])
  384. def test_enzo_example_c_with_degeneracy(self):
  385. # rescued from https://github.com/scipy/scipy/pull/218
  386. m = 20
  387. c = -np.ones(m)
  388. tmp = 2 * np.pi * np.arange(1, m + 1) / (m + 1)
  389. A_eq = np.vstack((np.cos(tmp) - 1, np.sin(tmp)))
  390. b_eq = [0, 0]
  391. res = linprog(c=c, A_eq=A_eq, b_eq=b_eq,
  392. method=self.method, options=self.options)
  393. _assert_success(res, desired_fun=0, desired_x=np.zeros(m))
  394. def test_enzo_example_c_with_unboundedness(self):
  395. # rescued from https://github.com/scipy/scipy/pull/218
  396. m = 50
  397. c = -np.ones(m)
  398. tmp = 2 * np.pi * np.arange(m) / (m + 1)
  399. A_eq = np.vstack((np.cos(tmp) - 1, np.sin(tmp)))
  400. b_eq = [0, 0]
  401. res = linprog(c=c, A_eq=A_eq, b_eq=b_eq,
  402. method=self.method, options=self.options)
  403. _assert_unbounded(res)
  404. def test_enzo_example_c_with_infeasibility(self):
  405. # rescued from https://github.com/scipy/scipy/pull/218
  406. m = 50
  407. c = -np.ones(m)
  408. tmp = 2 * np.pi * np.arange(m) / (m + 1)
  409. A_eq = np.vstack((np.cos(tmp) - 1, np.sin(tmp)))
  410. b_eq = [1, 1]
  411. o = {key: self.options[key] for key in self.options}
  412. o["presolve"] = False
  413. res = linprog(c=c, A_eq=A_eq, b_eq=b_eq, method=self.method,
  414. options=o)
  415. _assert_infeasible(res)
  416. def test_unknown_options(self):
  417. c = np.array([-3, -2])
  418. A_ub = [[2, 1], [1, 1], [1, 0]]
  419. b_ub = [10, 8, 4]
  420. def f(c, A_ub=None, b_ub=None, A_eq=None, b_eq=None, bounds=None,
  421. options={}):
  422. linprog(c, A_ub, b_ub, A_eq, b_eq, bounds, method=self.method,
  423. options=options)
  424. o = {key: self.options[key] for key in self.options}
  425. o['spam'] = 42
  426. _assert_warns(OptimizeWarning, f,
  427. c, A_ub=A_ub, b_ub=b_ub, options=o)
  428. def test_no_constraints(self):
  429. res = linprog([-1, -2], method=self.method, options=self.options)
  430. _assert_unbounded(res)
  431. def test_simple_bounds(self):
  432. res = linprog([1, 2], bounds=(1, 2),
  433. method=self.method, options=self.options)
  434. _assert_success(res, desired_x=[1, 1])
  435. res = linprog([1, 2], bounds=[(1, 2), (1, 2)],
  436. method=self.method, options=self.options)
  437. _assert_success(res, desired_x=[1, 1])
  438. def test_invalid_inputs(self):
  439. def f(c, A_ub=None, b_ub=None, A_eq=None, b_eq=None, bounds=None):
  440. linprog(c, A_ub, b_ub, A_eq, b_eq, bounds,
  441. method=self.method, options=self.options)
  442. for bad_bound in [[(5, 0), (1, 2), (3, 4)],
  443. [(1, 2), (3, 4)],
  444. [(1, 2), (3, 4), (3, 4, 5)],
  445. [(1, 2), (np.inf, np.inf), (3, 4)],
  446. [(1, 2), (-np.inf, -np.inf), (3, 4)],
  447. ]:
  448. assert_raises(ValueError, f, [1, 2, 3], bounds=bad_bound)
  449. assert_raises(ValueError, f, [1, 2], A_ub=[[1, 2]], b_ub=[1, 2])
  450. assert_raises(ValueError, f, [1, 2], A_ub=[[1]], b_ub=[1])
  451. assert_raises(ValueError, f, [1, 2], A_eq=[[1, 2]], b_eq=[1, 2])
  452. assert_raises(ValueError, f, [1, 2], A_eq=[[1]], b_eq=[1])
  453. assert_raises(ValueError, f, [1, 2], A_eq=[1], b_eq=1)
  454. if ("_sparse_presolve" in self.options and
  455. self.options["_sparse_presolve"]):
  456. return
  457. # this test doesn't make sense for sparse presolve
  458. # there aren't 3D sparse matrices
  459. assert_raises(ValueError, f, [1, 2], A_ub=np.zeros((1, 1, 3)), b_eq=1)
  460. def test_basic_artificial_vars(self):
  461. # Test if linprog succeeds when at the end of Phase 1 some artificial
  462. # variables remain basic, and the row in T corresponding to the
  463. # artificial variables is not all zero.
  464. c = np.array([-0.1, -0.07, 0.004, 0.004, 0.004, 0.004])
  465. A_ub = np.array([[1.0, 0, 0, 0, 0, 0], [-1.0, 0, 0, 0, 0, 0],
  466. [0, -1.0, 0, 0, 0, 0], [0, 1.0, 0, 0, 0, 0],
  467. [1.0, 1.0, 0, 0, 0, 0]])
  468. b_ub = np.array([3.0, 3.0, 3.0, 3.0, 20.0])
  469. A_eq = np.array([[1.0, 0, -1, 1, -1, 1], [0, -1.0, -1, 1, -1, 1]])
  470. b_eq = np.array([0, 0])
  471. res = linprog(c, A_ub=A_ub, b_ub=b_ub, A_eq=A_eq, b_eq=b_eq,
  472. method=self.method, options=self.options)
  473. _assert_success(res, desired_fun=0, desired_x=np.zeros_like(c),
  474. atol=2e-6)
  475. def test_empty_constraint_2(self):
  476. res = linprog([1, -1, 1, -1],
  477. bounds=[(0, np.inf), (-np.inf, 0), (-1, 1), (-1, 1)],
  478. method=self.method, options=self.options)
  479. _assert_success(res, desired_x=[0, 0, -1, 1], desired_fun=-2)
  480. def test_zero_row_2(self):
  481. A_eq = [[0, 0, 0], [1, 1, 1], [0, 0, 0]]
  482. b_eq = [0, 3, 0]
  483. c = [1, 2, 3]
  484. res = linprog(c=c, A_eq=A_eq, b_eq=b_eq,
  485. method=self.method, options=self.options)
  486. _assert_success(res, desired_fun=3)
  487. def test_zero_row_4(self):
  488. A_ub = [[0, 0, 0], [1, 1, 1], [0, 0, 0]]
  489. b_ub = [0, 3, 0]
  490. c = [1, 2, 3]
  491. res = linprog(c=c, A_ub=A_ub, b_ub=b_ub,
  492. method=self.method, options=self.options)
  493. _assert_success(res, desired_fun=0)
  494. def test_zero_column_1(self):
  495. m, n = 3, 4
  496. np.random.seed(0)
  497. c = np.random.rand(n)
  498. c[1] = 1
  499. A_eq = np.random.rand(m, n)
  500. A_eq[:, 1] = 0
  501. b_eq = np.random.rand(m)
  502. A_ub = [[1, 0, 1, 1]]
  503. b_ub = 3
  504. res = linprog(c, A_ub, b_ub, A_eq, b_eq,
  505. bounds=[(-10, 10), (-10, 10),
  506. (-10, None), (None, None)],
  507. method=self.method, options=self.options)
  508. _assert_success(res, desired_fun=-9.7087836730413404)
  509. def test_singleton_row_eq_2(self):
  510. c = [1, 1, 1, 2]
  511. A_eq = [[1, 0, 0, 0], [0, 2, 0, 0], [1, 0, 0, 0], [1, 1, 1, 1]]
  512. b_eq = [1, 2, 1, 4]
  513. res = linprog(c, A_eq=A_eq, b_eq=b_eq,
  514. method=self.method, options=self.options)
  515. _assert_success(res, desired_fun=4)
  516. def test_singleton_row_ub_2(self):
  517. c = [1, 1, 1, 2]
  518. A_ub = [[1, 0, 0, 0], [0, 2, 0, 0], [-1, 0, 0, 0], [1, 1, 1, 1]]
  519. b_ub = [1, 2, -0.5, 4]
  520. res = linprog(c, A_ub=A_ub, b_ub=b_ub,
  521. bounds=[(None, None), (0, None), (0, None), (0, None)],
  522. method=self.method, options=self.options)
  523. _assert_success(res, desired_fun=0.5)
  524. def test_remove_redundancy_infeasibility(self):
  525. m, n = 10, 10
  526. c = np.random.rand(n)
  527. A0 = np.random.rand(m, n)
  528. b0 = np.random.rand(m)
  529. A0[-1, :] = 2 * A0[-2, :]
  530. b0[-1] *= -1
  531. with suppress_warnings() as sup:
  532. sup.filter(OptimizeWarning, "A_eq does not appear...")
  533. res = linprog(c, A_eq=A0, b_eq=b0,
  534. method=self.method, options=self.options)
  535. _assert_infeasible(res)
  536. def test_bounded_below_only(self):
  537. A = np.eye(3)
  538. b = np.array([1, 2, 3])
  539. c = np.ones(3)
  540. res = linprog(c, A_eq=A, b_eq=b, bounds=(0.5, np.inf),
  541. method=self.method, options=self.options)
  542. _assert_success(res, desired_x=b, desired_fun=np.sum(b))
  543. def test_bounded_above_only(self):
  544. A = np.eye(3)
  545. b = np.array([1, 2, 3])
  546. c = np.ones(3)
  547. res = linprog(c, A_eq=A, b_eq=b, bounds=(-np.inf, 4),
  548. method=self.method, options=self.options)
  549. _assert_success(res, desired_x=b, desired_fun=np.sum(b))
  550. def test_unbounded_below_and_above(self):
  551. A = np.eye(3)
  552. b = np.array([1, 2, 3])
  553. c = np.ones(3)
  554. res = linprog(c, A_eq=A, b_eq=b, bounds=(-np.inf, np.inf),
  555. method=self.method, options=self.options)
  556. _assert_success(res, desired_x=b, desired_fun=np.sum(b))
  557. def test_bounds_equal_but_infeasible(self):
  558. c = [-4, 1]
  559. A_ub = [[7, -2], [0, 1], [2, -2]]
  560. b_ub = [14, 0, 3]
  561. bounds = [(2, 2), (0, None)]
  562. res = linprog(c=c, A_ub=A_ub, b_ub=b_ub, bounds=bounds,
  563. method=self.method, options=self.options)
  564. _assert_infeasible(res)
  565. def test_bounds_equal_but_infeasible2(self):
  566. c = [-4, 1]
  567. A_eq = [[7, -2], [0, 1], [2, -2]]
  568. b_eq = [14, 0, 3]
  569. bounds = [(2, 2), (0, None)]
  570. res = linprog(c=c, A_eq=A_eq, b_eq=b_eq, bounds=bounds,
  571. method=self.method, options=self.options)
  572. _assert_infeasible(res)
  573. def test_empty_constraint_1(self):
  574. res = linprog([-1, 1, -1, 1],
  575. bounds=[(0, np.inf), (-np.inf, 0), (-1, 1), (-1, 1)],
  576. method=self.method, options=self.options)
  577. _assert_unbounded(res)
  578. # Infeasibility detected in presolve requiring no iterations
  579. # if presolve is not used nit > 0 is expected.
  580. n = 0 if self.options.get('presolve', True) else 2
  581. assert_equal(res.nit, n)
  582. def test_singleton_row_eq_1(self):
  583. c = [1, 1, 1, 2]
  584. A_eq = [[1, 0, 0, 0], [0, 2, 0, 0], [1, 0, 0, 0], [1, 1, 1, 1]]
  585. b_eq = [1, 2, 2, 4]
  586. res = linprog(c, A_eq=A_eq, b_eq=b_eq,
  587. method=self.method, options=self.options)
  588. _assert_infeasible(res)
  589. # Infeasibility detected in presolve requiring no iterations
  590. # if presolve is not used nit > 0 is expected.
  591. n = 0 if self.options.get('presolve', True) else 3
  592. assert_equal(res.nit, n)
  593. def test_singleton_row_ub_1(self):
  594. c = [1, 1, 1, 2]
  595. A_ub = [[1, 0, 0, 0], [0, 2, 0, 0], [-1, 0, 0, 0], [1, 1, 1, 1]]
  596. b_ub = [1, 2, -2, 4]
  597. res = linprog(c, A_ub=A_ub, b_ub=b_ub,
  598. bounds=[(None, None), (0, None), (0, None), (0, None)],
  599. method=self.method, options=self.options)
  600. _assert_infeasible(res)
  601. # Infeasibility detected in presolve requiring no iterations
  602. # if presolve is not used nit > 0 is expected.
  603. n = 0 if self.options.get('presolve', True) else 3
  604. assert_equal(res.nit, n)
  605. def test_zero_column_2(self):
  606. np.random.seed(0)
  607. m, n = 2, 4
  608. c = np.random.rand(n)
  609. c[1] = -1
  610. A_eq = np.random.rand(m, n)
  611. A_eq[:, 1] = 0
  612. b_eq = np.random.rand(m)
  613. A_ub = np.random.rand(m, n)
  614. A_ub[:, 1] = 0
  615. b_ub = np.random.rand(m)
  616. res = linprog(c, A_ub, b_ub, A_eq, b_eq, bounds=(None, None),
  617. method=self.method, options=self.options)
  618. _assert_unbounded(res)
  619. # Infeasibility detected in presolve requiring no iterations
  620. # if presolve is not used nit > 0 is expected.
  621. n = 0 if self.options.get('presolve', True) else 5
  622. assert_equal(res.nit, n)
  623. def test_zero_row_1(self):
  624. m, n = 2, 4
  625. c = np.random.rand(n)
  626. A_eq = np.random.rand(m, n)
  627. A_eq[0, :] = 0
  628. b_eq = np.random.rand(m)
  629. res = linprog(c=c, A_eq=A_eq, b_eq=b_eq,
  630. method=self.method, options=self.options)
  631. _assert_infeasible(res)
  632. # Infeasibility detected in presolve requiring no iterations
  633. # if presolve is not used nit > 0 is expected.
  634. n = 0 if self.options.get('presolve', True) else 1
  635. assert_equal(res.nit, n)
  636. def test_zero_row_3(self):
  637. # detected in presolve?
  638. m, n = 2, 4
  639. c = np.random.rand(n)
  640. A_ub = np.random.rand(m, n)
  641. A_ub[0, :] = 0
  642. b_ub = -np.random.rand(m)
  643. res = linprog(c=c, A_ub=A_ub, b_ub=b_ub,
  644. method=self.method, options=self.options)
  645. _assert_infeasible(res)
  646. assert_equal(res.nit, 0)
  647. def test_infeasible_ub(self):
  648. c = [1]
  649. A_ub = [[2]]
  650. b_ub = 4
  651. bounds = (5, 6)
  652. res = linprog(c=c, A_ub=A_ub, b_ub=b_ub, bounds=bounds,
  653. method=self.method, options=self.options)
  654. _assert_infeasible(res)
  655. # Infeasibility detected in presolve requiring no iterations
  656. # if presolve is not used nit > 0 is expected.
  657. n = 0 if self.options.get('presolve', True) else 1
  658. assert_equal(res.nit, n)
  659. def test_type_error(self):
  660. c = [1]
  661. A_eq = [[1]]
  662. b_eq = "hello"
  663. assert_raises(TypeError, linprog,
  664. c, A_eq=A_eq, b_eq=b_eq,
  665. method=self.method, options=self.options)
  666. def test_equal_bounds_no_presolve(self):
  667. # There was a bug when a lower and upper bound were equal but
  668. # presolve was not on to eliminate the variable. The bound
  669. # was being converted to an equality constraint, but the bound
  670. # was not eliminated, leading to issues in postprocessing.
  671. c = [1, 2]
  672. A_ub = [[1, 2], [1.1, 2.2]]
  673. b_ub = [4, 8]
  674. bounds = [(1, 2), (2, 2)]
  675. o = {key: self.options[key] for key in self.options}
  676. o["presolve"] = False
  677. res = linprog(c=c, A_ub=A_ub, b_ub=b_ub, bounds=bounds,
  678. method=self.method, options=o)
  679. _assert_infeasible(res)
  680. def test_unbounded_below_no_presolve_corrected(self):
  681. c = [1]
  682. bounds = [(None, 1)]
  683. o = {key: self.options[key] for key in self.options}
  684. o["presolve"] = False
  685. res = linprog(c=c, bounds=bounds,
  686. method=self.method,
  687. options=o)
  688. _assert_unbounded(res)
  689. def test_unbounded_no_nontrivial_constraints_1(self):
  690. """
  691. Test whether presolve pathway for detecting unboundedness after
  692. constraint elimination is working.
  693. """
  694. c = np.array([0, 0, 0, 1, -1, -1])
  695. A = np.array([[1, 0, 0, 0, 0, 0],
  696. [0, 1, 0, 0, 0, 0],
  697. [0, 0, 0, 0, 0, -1]])
  698. b = np.array([2, -2, 0])
  699. bounds = [(None, None), (None, None), (None, None),
  700. (-1, 1), (-1, 1), (0, None)]
  701. res = linprog(c, A, b, None, None, bounds, method=self.method,
  702. options=self.options)
  703. _assert_unbounded(res)
  704. assert_equal(res.x[-1], np.inf)
  705. assert_equal(res.message[:36], "The problem is (trivially) unbounded")
  706. def test_unbounded_no_nontrivial_constraints_2(self):
  707. """
  708. Test whether presolve pathway for detecting unboundedness after
  709. constraint elimination is working.
  710. """
  711. c = np.array([0, 0, 0, 1, -1, 1])
  712. A = np.array([[1, 0, 0, 0, 0, 0],
  713. [0, 1, 0, 0, 0, 0],
  714. [0, 0, 0, 0, 0, 1]])
  715. b = np.array([2, -2, 0])
  716. bounds = [(None, None), (None, None), (None, None),
  717. (-1, 1), (-1, 1), (None, 0)]
  718. res = linprog(c, A, b, None, None, bounds, method=self.method,
  719. options=self.options)
  720. _assert_unbounded(res)
  721. assert_equal(res.x[-1], -np.inf)
  722. assert_equal(res.message[:36], "The problem is (trivially) unbounded")
  723. def test_bug_5400(self):
  724. # https://github.com/scipy/scipy/issues/5400
  725. bounds = [
  726. (0, None),
  727. (0, 100), (0, 100), (0, 100), (0, 100), (0, 100), (0, 100),
  728. (0, 900), (0, 900), (0, 900), (0, 900), (0, 900), (0, 900),
  729. (0, None), (0, None), (0, None), (0, None), (0, None), (0, None)]
  730. f = 1 / 9
  731. g = -1e4
  732. h = -3.1
  733. A_ub = np.array([
  734. [1, -2.99, 0, 0, -3, 0, 0, 0, -1, -1, 0, -1, -1, 1, 1, 0, 0, 0, 0],
  735. [1, 0, -2.9, h, 0, -3, 0, -1, 0, 0, -1, 0, -1, 0, 0, 1, 1, 0, 0],
  736. [1, 0, 0, h, 0, 0, -3, -1, -1, 0, -1, -1, 0, 0, 0, 0, 0, 1, 1],
  737. [0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
  738. [0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
  739. [0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
  740. [0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
  741. [0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
  742. [0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
  743. [0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
  744. [0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
  745. [0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0],
  746. [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0],
  747. [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0],
  748. [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0],
  749. [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, -1, 0, 0, 0, 0, 0],
  750. [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, -1, 0, 0, 0, 0],
  751. [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, -1, 0, 0, 0],
  752. [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, -1, 0, 0],
  753. [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, -1, 0],
  754. [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, -1],
  755. [0, 1.99, -1, -1, 0, 0, 0, -1, f, f, 0, 0, 0, g, 0, 0, 0, 0, 0],
  756. [0, 0, 0, 0, 2, -1, -1, 0, 0, 0, -1, f, f, 0, g, 0, 0, 0, 0],
  757. [0, -1, 1.9, 2.1, 0, 0, 0, f, -1, -1, 0, 0, 0, 0, 0, g, 0, 0, 0],
  758. [0, 0, 0, 0, -1, 2, -1, 0, 0, 0, f, -1, f, 0, 0, 0, g, 0, 0],
  759. [0, -1, -1, 2.1, 0, 0, 0, f, f, -1, 0, 0, 0, 0, 0, 0, 0, g, 0],
  760. [0, 0, 0, 0, -1, -1, 2, 0, 0, 0, f, f, -1, 0, 0, 0, 0, 0, g]])
  761. b_ub = np.array([
  762. 0.0, 0, 0, 100, 100, 100, 100, 100, 100, 900, 900, 900, 900, 900,
  763. 900, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0])
  764. c = np.array([-1.0, 1, 1, 1, 1, 1, 1, 1, 1,
  765. 1, 1, 1, 1, 0, 0, 0, 0, 0, 0])
  766. if self.method == 'simplex':
  767. with pytest.warns(OptimizeWarning):
  768. res = linprog(c, A_ub, b_ub, bounds=bounds,
  769. method=self.method, options=self.options)
  770. else:
  771. res = linprog(c, A_ub, b_ub, bounds=bounds,
  772. method=self.method, options=self.options)
  773. _assert_success(res, desired_fun=-106.63507541835018)
  774. def test_issue_6139(self):
  775. # Linprog(method='simplex') fails to find a basic feasible solution
  776. # if phase 1 pseudo-objective function is outside the provided tol.
  777. # https://github.com/scipy/scipy/issues/6139
  778. # Note: This is not strictly a bug as the default tolerance determines
  779. # if a result is "close enough" to zero and should not be expected
  780. # to work for all cases.
  781. c = np.array([1, 1, 1])
  782. A_eq = np.array([[1., 0., 0.], [-1000., 0., - 1000.]])
  783. b_eq = np.array([5.00000000e+00, -1.00000000e+04])
  784. A_ub = -np.array([[0., 1000000., 1010000.]])
  785. b_ub = -np.array([10000000.])
  786. bounds = (None, None)
  787. res = linprog(
  788. c, A_ub, b_ub, A_eq, b_eq, method=self.method,
  789. bounds=bounds, options=self.options
  790. )
  791. _assert_success(
  792. res, desired_fun=14.95, desired_x=np.array([5, 4.95, 5])
  793. )
  794. def test_bug_6690(self):
  795. # SciPy violates bound constraint despite result status being success
  796. # when the simplex method is used.
  797. # https://github.com/scipy/scipy/issues/6690
  798. A_eq = np.array([[0, 0, 0, 0.93, 0, 0.65, 0, 0, 0.83, 0]])
  799. b_eq = np.array([0.9626])
  800. A_ub = np.array([
  801. [0, 0, 0, 1.18, 0, 0, 0, -0.2, 0, -0.22],
  802. [0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
  803. [0, 0, 0, 0.43, 0, 0, 0, 0, 0, 0],
  804. [0, -1.22, -0.25, 0, 0, 0, -2.06, 0, 0, 1.37],
  805. [0, 0, 0, 0, 0, 0, 0, -0.25, 0, 0]
  806. ])
  807. b_ub = np.array([0.615, 0, 0.172, -0.869, -0.022])
  808. bounds = np.array([
  809. [-0.84, -0.97, 0.34, 0.4, -0.33, -0.74, 0.47, 0.09, -1.45, -0.73],
  810. [0.37, 0.02, 2.86, 0.86, 1.18, 0.5, 1.76, 0.17, 0.32, -0.15]
  811. ]).T
  812. c = np.array([
  813. -1.64, 0.7, 1.8, -1.06, -1.16, 0.26, 2.13, 1.53, 0.66, 0.28
  814. ])
  815. with suppress_warnings() as sup:
  816. sup.filter(OptimizeWarning, "Solving system with option 'sym_pos'")
  817. res = linprog(
  818. c, A_ub=A_ub, b_ub=b_ub, A_eq=A_eq, b_eq=b_eq,
  819. bounds=bounds, method=self.method, options=self.options
  820. )
  821. desired_fun = -1.19099999999
  822. desired_x = np.array([
  823. 0.3700, -0.9700, 0.3400, 0.4000, 1.1800,
  824. 0.5000, 0.4700, 0.0900, 0.3200, -0.7300
  825. ])
  826. _assert_success(
  827. res,
  828. desired_fun=desired_fun,
  829. desired_x=desired_x
  830. )
  831. # Add small tol value to ensure arrays are less than or equal.
  832. atol = 1e-6
  833. assert_array_less(bounds[:, 0] - atol, res.x)
  834. assert_array_less(res.x, bounds[:, 1] + atol)
  835. def test_bug_7044(self):
  836. # linprog fails to identify correct constraints with simplex method
  837. # leading to a non-optimal solution if A is rank-deficient.
  838. # https://github.com/scipy/scipy/issues/7044
  839. A, b, c, N = magic_square(3)
  840. with suppress_warnings() as sup:
  841. sup.filter(OptimizeWarning, "A_eq does not appear...")
  842. res = linprog(c, A_eq=A, b_eq=b,
  843. method=self.method, options=self.options)
  844. desired_fun = 1.730550597
  845. _assert_success(res, desired_fun=desired_fun)
  846. assert_allclose(A.dot(res.x), b)
  847. assert_array_less(np.zeros(res.x.size) - 1e-5, res.x)
  848. def test_issue_7237(self):
  849. # https://github.com/scipy/scipy/issues/7237
  850. # The simplex method sometimes "explodes" if the pivot value is very
  851. # close to zero.
  852. c = np.array([-1, 0, 0, 0, 0, 0, 0, 0, 0])
  853. A_ub = np.array([
  854. [1., -724., 911., -551., -555., -896., 478., -80., -293.],
  855. [1., 566., 42., 937.,233., 883., 392., -909., 57.],
  856. [1., -208., -894., 539., 321., 532., -924., 942., 55.],
  857. [1., 857., -859., 83., 462., -265., -971., 826., 482.],
  858. [1., 314., -424., 245., -424., 194., -443., -104., -429.],
  859. [1., 540., 679., 361., 149., -827., 876., 633., 302.],
  860. [0., -1., -0., -0., -0., -0., -0., -0., -0.],
  861. [0., -0., -1., -0., -0., -0., -0., -0., -0.],
  862. [0., -0., -0., -1., -0., -0., -0., -0., -0.],
  863. [0., -0., -0., -0., -1., -0., -0., -0., -0.],
  864. [0., -0., -0., -0., -0., -1., -0., -0., -0.],
  865. [0., -0., -0., -0., -0., -0., -1., -0., -0.],
  866. [0., -0., -0., -0., -0., -0., -0., -1., -0.],
  867. [0., -0., -0., -0., -0., -0., -0., -0., -1.],
  868. [0., 1., 0., 0., 0., 0., 0., 0., 0.],
  869. [0., 0., 1., 0., 0., 0., 0., 0., 0.],
  870. [0., 0., 0., 1., 0., 0., 0., 0., 0.],
  871. [0., 0., 0., 0., 1., 0., 0., 0., 0.],
  872. [0., 0., 0., 0., 0., 1., 0., 0., 0.],
  873. [0., 0., 0., 0., 0., 0., 1., 0., 0.],
  874. [0., 0., 0., 0., 0., 0., 0., 1., 0.],
  875. [0., 0., 0., 0., 0., 0., 0., 0., 1.]
  876. ])
  877. b_ub = np.array([
  878. 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.,
  879. 0., 0., 0., 1., 1., 1., 1., 1., 1., 1., 1.])
  880. A_eq = np.array([[0., 1., 1., 1., 1., 1., 1., 1., 1.]])
  881. b_eq = np.array([[1.]])
  882. bounds = [(None, None)] * 9
  883. res = linprog(c, A_ub=A_ub, b_ub=b_ub, A_eq=A_eq, b_eq=b_eq,
  884. bounds=bounds, method=self.method, options=self.options)
  885. _assert_success(res, desired_fun=108.568535, atol=1e-6)
  886. def test_issue_8174(self):
  887. # https://github.com/scipy/scipy/issues/8174
  888. # The simplex method sometimes "explodes" if the pivot value is very
  889. # close to zero.
  890. A_ub = np.array([
  891. [22714, 1008, 13380, -2713.5, -1116],
  892. [-4986, -1092, -31220, 17386.5, 684],
  893. [-4986, 0, 0, -2713.5, 0],
  894. [22714, 0, 0, 17386.5, 0]])
  895. b_ub = np.zeros(A_ub.shape[0])
  896. c = -np.ones(A_ub.shape[1])
  897. bounds = [(0,1)] * A_ub.shape[1]
  898. res = linprog(c=c, A_ub=A_ub, b_ub=b_ub, bounds=bounds,
  899. options=self.options, method=self.method)
  900. def test_issue_8174_stackoverflow(self):
  901. # Test supplementary example from issue 8174.
  902. # https://github.com/scipy/scipy/issues/8174
  903. # https://stackoverflow.com/questions/47717012/linprog-in-scipy-optimize-checking-solution
  904. c = np.array([1, 0, 0, 0, 0, 0, 0])
  905. A_ub = -np.identity(7)
  906. b_ub = np.array([[-2],[-2],[-2],[-2],[-2],[-2],[-2]])
  907. A_eq = np.array([
  908. [1, 1, 1, 1, 1, 1, 0],
  909. [0.3, 1.3, 0.9, 0, 0, 0, -1],
  910. [0.3, 0, 0, 0, 0, 0, -2/3],
  911. [0, 0.65, 0, 0, 0, 0, -1/15],
  912. [0, 0, 0.3, 0, 0, 0, -1/15]
  913. ])
  914. b_eq = np.array([[100],[0],[0],[0],[0]])
  915. with pytest.warns(OptimizeWarning):
  916. res = linprog(
  917. c=c, A_ub=A_ub, b_ub=b_ub, A_eq=A_eq, b_eq=b_eq,
  918. method=self.method, options=self.options
  919. )
  920. _assert_success(res, desired_fun=43.3333333331385)
  921. def test_bug_8662(self):
  922. # scipy.linprog returns incorrect optimal result for constraints using
  923. # default bounds, but, correct if boundary condition as constraint.
  924. # https://github.com/scipy/scipy/issues/8662
  925. c = [-10, 10, 6, 3]
  926. A = [
  927. [8, -8, -4, 6],
  928. [-8, 8, 4, -6],
  929. [-4, 4, 8, -4],
  930. [3, -3, -3, -10]
  931. ]
  932. b = [9, -9, -9, -4]
  933. bounds = [(0, None), (0, None), (0, None), (0, None)]
  934. desired_fun = 36.0000000000
  935. res1 = linprog(c, A, b, bounds=bounds,
  936. method=self.method, options=self.options)
  937. # Set boundary condition as a constraint
  938. A.append([0, 0, -1, 0])
  939. b.append(0)
  940. bounds[2] = (None, None)
  941. res2 = linprog(c, A, b, bounds=bounds, method=self.method,
  942. options=self.options)
  943. rtol = 1e-5
  944. _assert_success(res1, desired_fun=desired_fun, rtol=rtol)
  945. _assert_success(res2, desired_fun=desired_fun, rtol=rtol)
  946. def test_bug_8663(self):
  947. A = [[0, -7]]
  948. b = [-6]
  949. c = [1, 5]
  950. bounds = [(0, None), (None, None)]
  951. res = linprog(c, A_eq=A, b_eq=b, bounds=bounds,
  952. method=self.method, options=self.options)
  953. _assert_success(res,
  954. desired_x=[0, 6./7],
  955. desired_fun=5*6./7)
  956. def test_bug_8664(self):
  957. # Weak test. Ideally should _detect infeasibility_ for all options.
  958. c = [4]
  959. A_ub = [[2], [5]]
  960. b_ub = [4, 4]
  961. A_eq = [[0], [-8], [9]]
  962. b_eq = [3, 2, 10]
  963. with suppress_warnings() as sup:
  964. sup.filter(RuntimeWarning)
  965. sup.filter(OptimizeWarning, "Solving system with option...")
  966. o = {key: self.options[key] for key in self.options}
  967. o["presolve"] = False
  968. res = linprog(c, A_ub, b_ub, A_eq, b_eq, options=o,
  969. method=self.method)
  970. assert_(not res.success, "incorrectly reported success")
  971. def test_bug_8973(self):
  972. """
  973. Test whether bug described at:
  974. https://github.com/scipy/scipy/issues/8973
  975. was fixed.
  976. """
  977. c = np.array([0, 0, 0, 1, -1])
  978. A = np.array([[1, 0, 0, 0, 0], [0, 1, 0, 0, 0]])
  979. b = np.array([2, -2])
  980. bounds = [(None, None), (None, None), (None, None), (-1, 1), (-1, 1)]
  981. res = linprog(c, A, b, None, None, bounds, method=self.method,
  982. options=self.options)
  983. _assert_success(res,
  984. desired_x=[2, -2, 0, -1, 1],
  985. desired_fun=-2)
  986. def test_bug_8973_2(self):
  987. """
  988. Additional test for:
  989. https://github.com/scipy/scipy/issues/8973
  990. suggested in
  991. https://github.com/scipy/scipy/pull/8985
  992. review by @antonior92
  993. """
  994. c = np.zeros(1)
  995. A = np.array([[1]])
  996. b = np.array([-2])
  997. bounds = (None, None)
  998. res = linprog(c, A, b, None, None, bounds, method=self.method,
  999. options=self.options)
  1000. _assert_success(res) # would not pass if solution is infeasible
  1001. class BaseTestLinprogSimplex(LinprogCommonTests):
  1002. method = "simplex"
  1003. class TestLinprogSimplexCommon(BaseTestLinprogSimplex):
  1004. options = {}
  1005. def test_callback(self):
  1006. # Check that callback is as advertised
  1007. last_cb = {}
  1008. def cb(res):
  1009. message = res.pop('message')
  1010. complete = res.pop('complete')
  1011. assert_(res.pop('phase') in (1, 2))
  1012. assert_(res.pop('status') in range(4))
  1013. assert_(isinstance(res.pop('nit'), int))
  1014. assert_(isinstance(complete, bool))
  1015. assert_(isinstance(message, str))
  1016. if complete:
  1017. last_cb['x'] = res['x']
  1018. last_cb['fun'] = res['fun']
  1019. last_cb['slack'] = res['slack']
  1020. last_cb['con'] = res['con']
  1021. c = np.array([-3, -2])
  1022. A_ub = [[2, 1], [1, 1], [1, 0]]
  1023. b_ub = [10, 8, 4]
  1024. res = linprog(c, A_ub=A_ub, b_ub=b_ub, callback=cb, method=self.method)
  1025. _assert_success(res, desired_fun=-18.0, desired_x=[2, 6])
  1026. assert_allclose(last_cb['fun'], res['fun'])
  1027. assert_allclose(last_cb['x'], res['x'])
  1028. assert_allclose(last_cb['con'], res['con'])
  1029. assert_allclose(last_cb['slack'], res['slack'])
  1030. def test_issue_7237(self):
  1031. with pytest.raises(ValueError):
  1032. super(TestLinprogSimplexCommon, self).test_issue_7237()
  1033. def test_issue_8174(self):
  1034. with pytest.warns(OptimizeWarning):
  1035. super(TestLinprogSimplexCommon, self).test_issue_8174()
  1036. class TestLinprogSimplexBland(BaseTestLinprogSimplex):
  1037. options = {'bland': True}
  1038. def test_bug_5400(self):
  1039. with pytest.raises(ValueError):
  1040. super(TestLinprogSimplexBland, self).test_bug_5400()
  1041. def test_issue_8174(self):
  1042. with pytest.warns(OptimizeWarning):
  1043. super(TestLinprogSimplexBland, self).test_issue_8174()
  1044. class TestLinprogSimplexNoPresolve(BaseTestLinprogSimplex):
  1045. options = {'presolve': False}
  1046. def test_issue_6139(self):
  1047. # Linprog(method='simplex') fails to find a basic feasible solution
  1048. # if phase 1 pseudo-objective function is outside the provided tol.
  1049. # https://github.com/scipy/scipy/issues/6139
  1050. # Without ``presolve`` eliminating such rows the result is incorrect.
  1051. with pytest.raises(ValueError):
  1052. return super(TestLinprogSimplexNoPresolve, self).test_issue_6139()
  1053. def test_issue_7237(self):
  1054. with pytest.raises(ValueError):
  1055. super(TestLinprogSimplexNoPresolve, self).test_issue_7237()
  1056. def test_issue_8174(self):
  1057. with pytest.warns(OptimizeWarning):
  1058. super(TestLinprogSimplexNoPresolve, self).test_issue_8174()
  1059. def test_issue_8174_stackoverflow(self):
  1060. # Test expects linprog to raise a warning during presolve.
  1061. # As ``'presolve'=False`` no warning should be raised.
  1062. # Despite not presolving the result is still correct.
  1063. with pytest.warns(OptimizeWarning) as redundant_warning:
  1064. super(TestLinprogSimplexNoPresolve, self).test_issue_8174()
  1065. def test_unbounded_no_nontrivial_constraints_1(self):
  1066. pytest.skip("Tests behavior specific to presolve")
  1067. def test_unbounded_no_nontrivial_constraints_2(self):
  1068. pytest.skip("Tests behavior specific to presolve")
  1069. class BaseTestLinprogIP(LinprogCommonTests):
  1070. method = "interior-point"
  1071. class TestLinprogIPSpecific(object):
  1072. method = "interior-point"
  1073. # the following tests don't need to be performed separately for
  1074. # sparse presolve, sparse after presolve, and dense
  1075. def test_unbounded_below_no_presolve_original(self):
  1076. # formerly caused segfault in TravisCI w/ "cholesky":True
  1077. c = [-1]
  1078. bounds = [(None, 1)]
  1079. res = linprog(c=c, bounds=bounds,
  1080. method=self.method,
  1081. options={"presolve": False, "cholesky": True})
  1082. _assert_success(res, desired_fun=-1)
  1083. def test_cholesky(self):
  1084. # Test with a rather large problem (400 variables,
  1085. # 40 constraints) generated by https://gist.github.com/denis-bz/8647461
  1086. # use cholesky factorization and triangular solves
  1087. A, b, c = lpgen_2d(20, 20)
  1088. res = linprog(c, A_ub=A, b_ub=b, method=self.method,
  1089. options={"cholesky": True}) # only for dense
  1090. _assert_success(res, desired_fun=-64.049494229)
  1091. def test_alternate_initial_point(self):
  1092. # Test with a rather large problem (400 variables,
  1093. # 40 constraints) generated by https://gist.github.com/denis-bz/8647461
  1094. # use "improved" initial point
  1095. A, b, c = lpgen_2d(20, 20)
  1096. with suppress_warnings() as sup:
  1097. sup.filter(RuntimeWarning, "scipy.linalg.solve\nIll...")
  1098. sup.filter(OptimizeWarning, "Solving system with option...")
  1099. res = linprog(c, A_ub=A, b_ub=b, method=self.method,
  1100. options={"ip": True, "disp": True})
  1101. # ip code is independent of sparse/dense
  1102. _assert_success(res, desired_fun=-64.049494229)
  1103. def test_maxiter(self):
  1104. # Test with a rather large problem (400 variables,
  1105. # 40 constraints) generated by https://gist.github.com/denis-bz/8647461
  1106. # test iteration limit
  1107. A, b, c = lpgen_2d(20, 20)
  1108. maxiter = np.random.randint(6) + 1 # problem takes 7 iterations
  1109. res = linprog(c, A_ub=A, b_ub=b, method=self.method,
  1110. options={"maxiter": maxiter})
  1111. # maxiter is independent of sparse/dense
  1112. assert_equal(res.status, 1)
  1113. assert_equal(res.nit, maxiter)
  1114. def test_disp(self):
  1115. # Test with a rather large problem (400 variables,
  1116. # 40 constraints) generated by https://gist.github.com/denis-bz/8647461
  1117. # test that display option does not break anything.
  1118. A, b, c = lpgen_2d(20, 20)
  1119. res = linprog(c, A_ub=A, b_ub=b, method=self.method,
  1120. options={"disp": True})
  1121. # disp is independent of sparse/dense
  1122. _assert_success(res, desired_fun=-64.049494229)
  1123. def test_callback(self):
  1124. def f():
  1125. pass
  1126. A = [[0, -7]]
  1127. b = [-6]
  1128. c = [1, 5]
  1129. bounds = [(0, None), (None, None)]
  1130. # Linprog should solve in presolve. As the interior-point method is
  1131. # not used the the callback should never be needed and no error
  1132. # returned
  1133. res = linprog(c, A_eq=A, b_eq=b, bounds=bounds,
  1134. method=self.method, callback=f)
  1135. _assert_success(res, desired_x=[0, 6./7], desired_fun=5*6./7)
  1136. # Without presolve the solver reverts to the interior-point method
  1137. # Interior-point currently does not implement callback functions.
  1138. with pytest.raises(NotImplementedError):
  1139. res = linprog(c, A_eq=A, b_eq=b, bounds=bounds, method=self.method,
  1140. callback=f, options={'presolve': False})
  1141. class TestLinprogIPSparse(BaseTestLinprogIP):
  1142. options = {"sparse": True}
  1143. @pytest.mark.xfail(reason='Fails with ATLAS, see gh-7877')
  1144. def test_bug_6690(self):
  1145. # Test defined in base class, but can't mark as xfail there
  1146. super(TestLinprogIPSparse, self).test_bug_6690()
  1147. def test_magic_square_sparse_no_presolve(self):
  1148. # test linprog with a problem with a rank-deficient A_eq matrix
  1149. A, b, c, N = magic_square(3)
  1150. with suppress_warnings() as sup:
  1151. sup.filter(MatrixRankWarning, "Matrix is exactly singular")
  1152. sup.filter(OptimizeWarning, "Solving system with option...")
  1153. o = {key: self.options[key] for key in self.options}
  1154. o["presolve"] = False
  1155. res = linprog(c, A_eq=A, b_eq=b, bounds=(0, 1),
  1156. options=o, method=self.method)
  1157. _assert_success(res, desired_fun=1.730550597)
  1158. def test_sparse_solve_options(self):
  1159. A, b, c, N = magic_square(3)
  1160. with suppress_warnings() as sup:
  1161. sup.filter(OptimizeWarning, "A_eq does not appear...")
  1162. sup.filter(OptimizeWarning, "Invalid permc_spec option")
  1163. o = {key: self.options[key] for key in self.options}
  1164. permc_specs = ('NATURAL', 'MMD_ATA', 'MMD_AT_PLUS_A',
  1165. 'COLAMD', 'ekki-ekki-ekki')
  1166. for permc_spec in permc_specs:
  1167. o["permc_spec"] = permc_spec
  1168. res = linprog(c, A_eq=A, b_eq=b, bounds=(0, 1),
  1169. method=self.method, options=o)
  1170. _assert_success(res, desired_fun=1.730550597)
  1171. class TestLinprogIPDense(BaseTestLinprogIP):
  1172. options = {"sparse": False}
  1173. class TestLinprogIPSparsePresolve(BaseTestLinprogIP):
  1174. options = {"sparse": True, "_sparse_presolve": True}
  1175. def test_enzo_example_c_with_infeasibility(self):
  1176. pytest.skip('_sparse_presolve=True incompatible with presolve=False')
  1177. @pytest.mark.xfail(reason='Fails with ATLAS, see gh-7877')
  1178. def test_bug_6690(self):
  1179. # Test defined in base class, but can't mark as xfail there
  1180. super(TestLinprogIPSparsePresolve, self).test_bug_6690()
  1181. def test_unknown_solver():
  1182. c = np.array([-3, -2])
  1183. A_ub = [[2, 1], [1, 1], [1, 0]]
  1184. b_ub = [10, 8, 4]
  1185. assert_raises(ValueError, linprog,
  1186. c, A_ub=A_ub, b_ub=b_ub, method='ekki-ekki-ekki')