Coverage for /usr/lib/python3/dist-packages/scipy/optimize/_cobyla_py.py: 19%

98 statements  

« prev     ^ index     » next       coverage.py v7.9.1, created at 2025-06-14 15:55 +0200

1""" 

2Interface to Constrained Optimization By Linear Approximation 

3 

4Functions 

5--------- 

6.. autosummary:: 

7 :toctree: generated/ 

8 

9 fmin_cobyla 

10 

11""" 

12 

13import functools 

14from threading import RLock 

15 

16import numpy as np 

17from scipy.optimize import _cobyla as cobyla 

18from ._optimize import (OptimizeResult, _check_unknown_options, 

19 _prepare_scalar_function) 

20try: 

21 from itertools import izip 

22except ImportError: 

23 izip = zip 

24 

25__all__ = ['fmin_cobyla'] 

26 

27# Workarund as _cobyla.minimize is not threadsafe 

28# due to an unknown f2py bug and can segfault, 

29# see gh-9658. 

30_module_lock = RLock() 

31def synchronized(func): 

32 @functools.wraps(func) 

33 def wrapper(*args, **kwargs): 

34 with _module_lock: 

35 return func(*args, **kwargs) 

36 return wrapper 

37 

38@synchronized 

39def fmin_cobyla(func, x0, cons, args=(), consargs=None, rhobeg=1.0, 

40 rhoend=1e-4, maxfun=1000, disp=None, catol=2e-4, 

41 *, callback=None): 

42 """ 

43 Minimize a function using the Constrained Optimization By Linear 

44 Approximation (COBYLA) method. This method wraps a FORTRAN 

45 implementation of the algorithm. 

46 

47 Parameters 

48 ---------- 

49 func : callable 

50 Function to minimize. In the form func(x, \\*args). 

51 x0 : ndarray 

52 Initial guess. 

53 cons : sequence 

54 Constraint functions; must all be ``>=0`` (a single function 

55 if only 1 constraint). Each function takes the parameters `x` 

56 as its first argument, and it can return either a single number or 

57 an array or list of numbers. 

58 args : tuple, optional 

59 Extra arguments to pass to function. 

60 consargs : tuple, optional 

61 Extra arguments to pass to constraint functions (default of None means 

62 use same extra arguments as those passed to func). 

63 Use ``()`` for no extra arguments. 

64 rhobeg : float, optional 

65 Reasonable initial changes to the variables. 

66 rhoend : float, optional 

67 Final accuracy in the optimization (not precisely guaranteed). This 

68 is a lower bound on the size of the trust region. 

69 disp : {0, 1, 2, 3}, optional 

70 Controls the frequency of output; 0 implies no output. 

71 maxfun : int, optional 

72 Maximum number of function evaluations. 

73 catol : float, optional 

74 Absolute tolerance for constraint violations. 

75 callback : callable, optional 

76 Called after each iteration, as ``callback(x)``, where ``x`` is the 

77 current parameter vector. 

78 

79 Returns 

80 ------- 

81 x : ndarray 

82 The argument that minimises `f`. 

83 

84 See also 

85 -------- 

86 minimize: Interface to minimization algorithms for multivariate 

87 functions. See the 'COBYLA' `method` in particular. 

88 

89 Notes 

90 ----- 

91 This algorithm is based on linear approximations to the objective 

92 function and each constraint. We briefly describe the algorithm. 

93 

94 Suppose the function is being minimized over k variables. At the 

95 jth iteration the algorithm has k+1 points v_1, ..., v_(k+1), 

96 an approximate solution x_j, and a radius RHO_j. 

97 (i.e., linear plus a constant) approximations to the objective 

98 function and constraint functions such that their function values 

99 agree with the linear approximation on the k+1 points v_1,.., v_(k+1). 

100 This gives a linear program to solve (where the linear approximations 

101 of the constraint functions are constrained to be non-negative). 

102 

103 However, the linear approximations are likely only good 

104 approximations near the current simplex, so the linear program is 

105 given the further requirement that the solution, which 

106 will become x_(j+1), must be within RHO_j from x_j. RHO_j only 

107 decreases, never increases. The initial RHO_j is rhobeg and the 

108 final RHO_j is rhoend. In this way COBYLA's iterations behave 

109 like a trust region algorithm. 

110 

111 Additionally, the linear program may be inconsistent, or the 

112 approximation may give poor improvement. For details about 

113 how these issues are resolved, as well as how the points v_i are 

114 updated, refer to the source code or the references below. 

115 

116 

117 References 

118 ---------- 

119 Powell M.J.D. (1994), "A direct search optimization method that models 

120 the objective and constraint functions by linear interpolation.", in 

121 Advances in Optimization and Numerical Analysis, eds. S. Gomez and 

122 J-P Hennart, Kluwer Academic (Dordrecht), pp. 51-67 

123 

124 Powell M.J.D. (1998), "Direct search algorithms for optimization 

125 calculations", Acta Numerica 7, 287-336 

126 

127 Powell M.J.D. (2007), "A view of algorithms for optimization without 

128 derivatives", Cambridge University Technical Report DAMTP 2007/NA03 

129 

130 

131 Examples 

132 -------- 

133 Minimize the objective function f(x,y) = x*y subject 

134 to the constraints x**2 + y**2 < 1 and y > 0:: 

135 

136 >>> def objective(x): 

137 ... return x[0]*x[1] 

138 ... 

139 >>> def constr1(x): 

140 ... return 1 - (x[0]**2 + x[1]**2) 

141 ... 

142 >>> def constr2(x): 

143 ... return x[1] 

144 ... 

145 >>> from scipy.optimize import fmin_cobyla 

146 >>> fmin_cobyla(objective, [0.0, 0.1], [constr1, constr2], rhoend=1e-7) 

147 array([-0.70710685, 0.70710671]) 

148 

149 The exact solution is (-sqrt(2)/2, sqrt(2)/2). 

150 

151 

152 

153 """ 

154 err = "cons must be a sequence of callable functions or a single"\ 

155 " callable function." 

156 try: 

157 len(cons) 

158 except TypeError as e: 

159 if callable(cons): 

160 cons = [cons] 

161 else: 

162 raise TypeError(err) from e 

163 else: 

164 for thisfunc in cons: 

165 if not callable(thisfunc): 

166 raise TypeError(err) 

167 

168 if consargs is None: 

169 consargs = args 

170 

171 # build constraints 

172 con = tuple({'type': 'ineq', 'fun': c, 'args': consargs} for c in cons) 

173 

174 # options 

175 opts = {'rhobeg': rhobeg, 

176 'tol': rhoend, 

177 'disp': disp, 

178 'maxiter': maxfun, 

179 'catol': catol, 

180 'callback': callback} 

181 

182 sol = _minimize_cobyla(func, x0, args, constraints=con, 

183 **opts) 

184 if disp and not sol['success']: 

185 print(f"COBYLA failed to find a solution: {sol.message}") 

186 return sol['x'] 

187 

188 

189@synchronized 

190def _minimize_cobyla(fun, x0, args=(), constraints=(), 

191 rhobeg=1.0, tol=1e-4, maxiter=1000, 

192 disp=False, catol=2e-4, callback=None, bounds=None, 

193 **unknown_options): 

194 """ 

195 Minimize a scalar function of one or more variables using the 

196 Constrained Optimization BY Linear Approximation (COBYLA) algorithm. 

197 

198 Options 

199 ------- 

200 rhobeg : float 

201 Reasonable initial changes to the variables. 

202 tol : float 

203 Final accuracy in the optimization (not precisely guaranteed). 

204 This is a lower bound on the size of the trust region. 

205 disp : bool 

206 Set to True to print convergence messages. If False, 

207 `verbosity` is ignored as set to 0. 

208 maxiter : int 

209 Maximum number of function evaluations. 

210 catol : float 

211 Tolerance (absolute) for constraint violations 

212 

213 """ 

214 _check_unknown_options(unknown_options) 

215 maxfun = maxiter 

216 rhoend = tol 

217 iprint = int(bool(disp)) 

218 

219 # check constraints 

220 if isinstance(constraints, dict): 

221 constraints = (constraints, ) 

222 

223 if bounds: 

224 i_lb = np.isfinite(bounds.lb) 

225 if np.any(i_lb): 

226 def lb_constraint(x, *args, **kwargs): 

227 return x[i_lb] - bounds.lb[i_lb] 

228 

229 constraints.append({'type': 'ineq', 'fun': lb_constraint}) 

230 

231 i_ub = np.isfinite(bounds.ub) 

232 if np.any(i_ub): 

233 def ub_constraint(x): 

234 return bounds.ub[i_ub] - x[i_ub] 

235 

236 constraints.append({'type': 'ineq', 'fun': ub_constraint}) 

237 

238 for ic, con in enumerate(constraints): 

239 # check type 

240 try: 

241 ctype = con['type'].lower() 

242 except KeyError as e: 

243 raise KeyError('Constraint %d has no type defined.' % ic) from e 

244 except TypeError as e: 

245 raise TypeError('Constraints must be defined using a ' 

246 'dictionary.') from e 

247 except AttributeError as e: 

248 raise TypeError("Constraint's type must be a string.") from e 

249 else: 

250 if ctype != 'ineq': 

251 raise ValueError("Constraints of type '%s' not handled by " 

252 "COBYLA." % con['type']) 

253 

254 # check function 

255 if 'fun' not in con: 

256 raise KeyError('Constraint %d has no function defined.' % ic) 

257 

258 # check extra arguments 

259 if 'args' not in con: 

260 con['args'] = () 

261 

262 # m is the total number of constraint values 

263 # it takes into account that some constraints may be vector-valued 

264 cons_lengths = [] 

265 for c in constraints: 

266 f = c['fun'](x0, *c['args']) 

267 try: 

268 cons_length = len(f) 

269 except TypeError: 

270 cons_length = 1 

271 cons_lengths.append(cons_length) 

272 m = sum(cons_lengths) 

273 

274 # create the ScalarFunction, cobyla doesn't require derivative function 

275 def _jac(x, *args): 

276 return None 

277 

278 sf = _prepare_scalar_function(fun, x0, args=args, jac=_jac) 

279 

280 def calcfc(x, con): 

281 f = sf.fun(x) 

282 i = 0 

283 for size, c in izip(cons_lengths, constraints): 

284 con[i: i + size] = c['fun'](x, *c['args']) 

285 i += size 

286 return f 

287 

288 def wrapped_callback(x): 

289 if callback is not None: 

290 callback(np.copy(x)) 

291 

292 info = np.zeros(4, np.float64) 

293 xopt, info = cobyla.minimize(calcfc, m=m, x=np.copy(x0), rhobeg=rhobeg, 

294 rhoend=rhoend, iprint=iprint, maxfun=maxfun, 

295 dinfo=info, callback=wrapped_callback) 

296 

297 if info[3] > catol: 

298 # Check constraint violation 

299 info[0] = 4 

300 

301 return OptimizeResult(x=xopt, 

302 status=int(info[0]), 

303 success=info[0] == 1, 

304 message={1: 'Optimization terminated successfully.', 

305 2: 'Maximum number of function evaluations ' 

306 'has been exceeded.', 

307 3: 'Rounding errors are becoming damaging ' 

308 'in COBYLA subroutine.', 

309 4: 'Did not converge to a solution ' 

310 'satisfying the constraints. See ' 

311 '`maxcv` for magnitude of violation.', 

312 5: 'NaN result encountered.' 

313 }.get(info[0], 'Unknown exit status.'), 

314 nfev=int(info[1]), 

315 fun=info[2], 

316 maxcv=info[3])