test_cobyla.py 6.7 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195
  1. import math
  2. import numpy as np
  3. from numpy.testing import assert_allclose, assert_array_almost_equal
  4. from scipy.optimize import (
  5. fmin_cobyla, minimize, Bounds, NonlinearConstraint, LinearConstraint,
  6. OptimizeResult
  7. )
  8. class TestCobyla:
  9. def setup_method(self):
  10. # The algorithm is very fragile on 32 bit, so unfortunately we need to start
  11. # very near the solution in order for the test to pass.
  12. self.x0 = [np.sqrt(25 - (2.0/3)**2), 2.0/3 + 1e-4]
  13. self.solution = [math.sqrt(25 - (2.0/3)**2), 2.0/3]
  14. self.opts = {'disp': 0, 'rhobeg': 1, 'tol': 1e-6,
  15. 'maxiter': 100}
  16. def fun(self, x):
  17. return x[0]**2 + abs(x[1])**3
  18. def con1(self, x):
  19. return x[0]**2 + x[1]**2 - 25
  20. def con2(self, x):
  21. return -self.con1(x)
  22. def test_simple(self):
  23. # use disp=True as smoke test for gh-8118
  24. x = fmin_cobyla(self.fun, self.x0, [self.con1, self.con2], rhobeg=1,
  25. rhoend=1e-5, maxfun=100, disp=1)
  26. assert_allclose(x, self.solution, atol=1e-4)
  27. def test_minimize_simple(self):
  28. class Callback:
  29. def __init__(self):
  30. self.n_calls = 0
  31. self.last_x = None
  32. def __call__(self, x):
  33. self.n_calls += 1
  34. self.last_x = x
  35. class CallbackNewSyntax:
  36. def __init__(self):
  37. self.n_calls = 0
  38. def __call__(self, intermediate_result):
  39. assert isinstance(intermediate_result, OptimizeResult)
  40. self.n_calls += 1
  41. callback = Callback()
  42. callback_new_syntax = CallbackNewSyntax()
  43. # Minimize with method='COBYLA'
  44. cons = (NonlinearConstraint(self.con1, 0, np.inf),
  45. {'type': 'ineq', 'fun': self.con2})
  46. sol = minimize(self.fun, self.x0, method='cobyla', constraints=cons,
  47. callback=callback, options=self.opts)
  48. sol_new = minimize(self.fun, self.x0, method='cobyla', constraints=cons,
  49. callback=callback_new_syntax, options=self.opts)
  50. assert_allclose(sol.x, self.solution, atol=1e-4)
  51. assert sol.success, sol.message
  52. assert sol.maxcv < 1e-5, sol
  53. assert sol.nfev < 70, sol
  54. assert sol.fun < self.fun(self.solution) + 1e-3, sol
  55. assert_array_almost_equal(
  56. sol.x,
  57. callback.last_x,
  58. decimal=5,
  59. err_msg="Last design vector sent to the callback is not equal to"
  60. " returned value.",
  61. )
  62. assert sol_new.success, sol_new.message
  63. assert sol.fun == sol_new.fun
  64. assert sol.maxcv == sol_new.maxcv
  65. assert sol.nfev == sol_new.nfev
  66. assert callback.n_calls == callback_new_syntax.n_calls, \
  67. "Callback is not called the same amount of times for old and new syntax."
  68. def test_minimize_constraint_violation(self):
  69. # We set up conflicting constraints so that the algorithm will be
  70. # guaranteed to end up with maxcv > 0.
  71. cons = ({'type': 'ineq', 'fun': lambda x: 4 - x},
  72. {'type': 'ineq', 'fun': lambda x: x - 5})
  73. sol = minimize(lambda x: x, [0], method='cobyla', constraints=cons,
  74. options={'catol': 0.6})
  75. assert sol.maxcv > 0.1
  76. assert sol.success
  77. sol = minimize(lambda x: x, [0], method='cobyla', constraints=cons,
  78. options={'catol': 0.4})
  79. assert sol.maxcv > 0.1
  80. assert not sol.success
  81. def test_f_target(self):
  82. f_target = 250
  83. sol = minimize(lambda x: x**2, [500], method='cobyla',
  84. options={'f_target': f_target})
  85. assert sol.status == 1
  86. assert sol.success
  87. assert sol.fun <= f_target
  88. def test_minimize_linear_constraints(self):
  89. constraints = LinearConstraint([1.0, 1.0], 1.0, 1.0)
  90. sol = minimize(
  91. self.fun,
  92. self.x0,
  93. method='cobyla',
  94. constraints=constraints,
  95. options=self.opts,
  96. )
  97. solution = [(4 - np.sqrt(7)) / 3, (np.sqrt(7) - 1) / 3]
  98. assert_allclose(sol.x, solution, atol=1e-4)
  99. assert sol.success, sol.message
  100. assert sol.maxcv < 1e-8, sol
  101. assert sol.nfev <= 100, sol
  102. assert sol.fun < self.fun(solution) + 1e-3, sol
  103. def test_vector_constraints():
  104. # test that fmin_cobyla and minimize can take a combination
  105. # of constraints, some returning a number and others an array
  106. def fun(x):
  107. return (x[0] - 1)**2 + (x[1] - 2.5)**2
  108. def fmin(x):
  109. return fun(x) - 1
  110. def cons1(x):
  111. a = np.array([[1, -2, 2], [-1, -2, 6], [-1, 2, 2]])
  112. return np.array([a[i, 0] * x[0] + a[i, 1] * x[1] +
  113. a[i, 2] for i in range(len(a))])
  114. def cons2(x):
  115. return x # identity, acts as bounds x > 0
  116. x0 = np.array([2, 0])
  117. cons_list = [fun, cons1, cons2]
  118. xsol = [1.4, 1.7]
  119. fsol = 0.8
  120. # testing fmin_cobyla
  121. sol = fmin_cobyla(fun, x0, cons_list, rhoend=1e-5)
  122. assert_allclose(sol, xsol, atol=1e-4)
  123. sol = fmin_cobyla(fun, x0, fmin, rhoend=1e-5)
  124. assert_allclose(fun(sol), 1, atol=1e-4)
  125. # testing minimize
  126. constraints = [{'type': 'ineq', 'fun': cons} for cons in cons_list]
  127. sol = minimize(fun, x0, constraints=constraints, tol=1e-5)
  128. assert_allclose(sol.x, xsol, atol=1e-4)
  129. assert sol.success, sol.message
  130. assert_allclose(sol.fun, fsol, atol=1e-4)
  131. constraints = {'type': 'ineq', 'fun': fmin}
  132. sol = minimize(fun, x0, constraints=constraints, tol=1e-5)
  133. assert_allclose(sol.fun, 1, atol=1e-4)
  134. class TestBounds:
  135. # Test cobyla support for bounds (only when used via `minimize`)
  136. # Invalid bounds is tested in
  137. # test_optimize.TestOptimizeSimple.test_minimize_invalid_bounds
  138. def test_basic(self):
  139. def f(x):
  140. return np.sum(x**2)
  141. lb = [-1, None, 1, None, -0.5]
  142. ub = [-0.5, -0.5, None, None, -0.5]
  143. bounds = [(a, b) for a, b in zip(lb, ub)]
  144. # these are converted to Bounds internally
  145. res = minimize(f, x0=[1, 2, 3, 4, 5], method='cobyla', bounds=bounds)
  146. ref = [-0.5, -0.5, 1, 0, -0.5]
  147. assert res.success
  148. assert_allclose(res.x, ref, atol=1e-3)
  149. def test_unbounded(self):
  150. def f(x):
  151. return np.sum(x**2)
  152. bounds = Bounds([-np.inf, -np.inf], [np.inf, np.inf])
  153. res = minimize(f, x0=[1, 2], method='cobyla', bounds=bounds)
  154. assert res.success
  155. assert_allclose(res.x, 0, atol=1e-3)
  156. bounds = Bounds([1, -np.inf], [np.inf, np.inf])
  157. res = minimize(f, x0=[1, 2], method='cobyla', bounds=bounds)
  158. assert res.success
  159. assert_allclose(res.x, [1, 0], atol=1e-3)