Coverage for /pythoncovmergedfiles/medio/medio/usr/local/lib/python3.8/site-packages/scipy/optimize/_differentiable_functions.py: 9%
364 statements
« prev ^ index » next coverage.py v7.3.2, created at 2023-12-12 06:31 +0000
1import numpy as np
2import scipy.sparse as sps
3from ._numdiff import approx_derivative, group_columns
4from ._hessian_update_strategy import HessianUpdateStrategy
5from scipy.sparse.linalg import LinearOperator
8FD_METHODS = ('2-point', '3-point', 'cs')
11class ScalarFunction:
12 """Scalar function and its derivatives.
14 This class defines a scalar function F: R^n->R and methods for
15 computing or approximating its first and second derivatives.
17 Parameters
18 ----------
19 fun : callable
20 evaluates the scalar function. Must be of the form ``fun(x, *args)``,
21 where ``x`` is the argument in the form of a 1-D array and ``args`` is
22 a tuple of any additional fixed parameters needed to completely specify
23 the function. Should return a scalar.
24 x0 : array-like
25 Provides an initial set of variables for evaluating fun. Array of real
26 elements of size (n,), where 'n' is the number of independent
27 variables.
28 args : tuple, optional
29 Any additional fixed parameters needed to completely specify the scalar
30 function.
31 grad : {callable, '2-point', '3-point', 'cs'}
32 Method for computing the gradient vector.
33 If it is a callable, it should be a function that returns the gradient
34 vector:
36 ``grad(x, *args) -> array_like, shape (n,)``
38 where ``x`` is an array with shape (n,) and ``args`` is a tuple with
39 the fixed parameters.
40 Alternatively, the keywords {'2-point', '3-point', 'cs'} can be used
41 to select a finite difference scheme for numerical estimation of the
42 gradient with a relative step size. These finite difference schemes
43 obey any specified `bounds`.
44 hess : {callable, '2-point', '3-point', 'cs', HessianUpdateStrategy}
45 Method for computing the Hessian matrix. If it is callable, it should
46 return the Hessian matrix:
48 ``hess(x, *args) -> {LinearOperator, spmatrix, array}, (n, n)``
50 where x is a (n,) ndarray and `args` is a tuple with the fixed
51 parameters. Alternatively, the keywords {'2-point', '3-point', 'cs'}
52 select a finite difference scheme for numerical estimation. Or, objects
53 implementing `HessianUpdateStrategy` interface can be used to
54 approximate the Hessian.
55 Whenever the gradient is estimated via finite-differences, the Hessian
56 cannot be estimated with options {'2-point', '3-point', 'cs'} and needs
57 to be estimated using one of the quasi-Newton strategies.
58 finite_diff_rel_step : None or array_like
59 Relative step size to use. The absolute step size is computed as
60 ``h = finite_diff_rel_step * sign(x0) * max(1, abs(x0))``, possibly
61 adjusted to fit into the bounds. For ``method='3-point'`` the sign
62 of `h` is ignored. If None then finite_diff_rel_step is selected
63 automatically,
64 finite_diff_bounds : tuple of array_like
65 Lower and upper bounds on independent variables. Defaults to no bounds,
66 (-np.inf, np.inf). Each bound must match the size of `x0` or be a
67 scalar, in the latter case the bound will be the same for all
68 variables. Use it to limit the range of function evaluation.
69 epsilon : None or array_like, optional
70 Absolute step size to use, possibly adjusted to fit into the bounds.
71 For ``method='3-point'`` the sign of `epsilon` is ignored. By default
72 relative steps are used, only if ``epsilon is not None`` are absolute
73 steps used.
75 Notes
76 -----
77 This class implements a memoization logic. There are methods `fun`,
78 `grad`, hess` and corresponding attributes `f`, `g` and `H`. The following
79 things should be considered:
81 1. Use only public methods `fun`, `grad` and `hess`.
82 2. After one of the methods is called, the corresponding attribute
83 will be set. However, a subsequent call with a different argument
84 of *any* of the methods may overwrite the attribute.
85 """
86 def __init__(self, fun, x0, args, grad, hess, finite_diff_rel_step,
87 finite_diff_bounds, epsilon=None):
88 if not callable(grad) and grad not in FD_METHODS:
89 raise ValueError(
90 f"`grad` must be either callable or one of {FD_METHODS}."
91 )
93 if not (callable(hess) or hess in FD_METHODS
94 or isinstance(hess, HessianUpdateStrategy)):
95 raise ValueError(
96 f"`hess` must be either callable, HessianUpdateStrategy"
97 f" or one of {FD_METHODS}."
98 )
100 if grad in FD_METHODS and hess in FD_METHODS:
101 raise ValueError("Whenever the gradient is estimated via "
102 "finite-differences, we require the Hessian "
103 "to be estimated using one of the "
104 "quasi-Newton strategies.")
106 # the astype call ensures that self.x is a copy of x0
107 self.x = np.atleast_1d(x0).astype(float)
108 self.n = self.x.size
109 self.nfev = 0
110 self.ngev = 0
111 self.nhev = 0
112 self.f_updated = False
113 self.g_updated = False
114 self.H_updated = False
116 self._lowest_x = None
117 self._lowest_f = np.inf
119 finite_diff_options = {}
120 if grad in FD_METHODS:
121 finite_diff_options["method"] = grad
122 finite_diff_options["rel_step"] = finite_diff_rel_step
123 finite_diff_options["abs_step"] = epsilon
124 finite_diff_options["bounds"] = finite_diff_bounds
125 if hess in FD_METHODS:
126 finite_diff_options["method"] = hess
127 finite_diff_options["rel_step"] = finite_diff_rel_step
128 finite_diff_options["abs_step"] = epsilon
129 finite_diff_options["as_linear_operator"] = True
131 # Function evaluation
132 def fun_wrapped(x):
133 self.nfev += 1
134 # Send a copy because the user may overwrite it.
135 # Overwriting results in undefined behaviour because
136 # fun(self.x) will change self.x, with the two no longer linked.
137 fx = fun(np.copy(x), *args)
138 # Make sure the function returns a true scalar
139 if not np.isscalar(fx):
140 try:
141 fx = np.asarray(fx).item()
142 except (TypeError, ValueError) as e:
143 raise ValueError(
144 "The user-provided objective function "
145 "must return a scalar value."
146 ) from e
148 if fx < self._lowest_f:
149 self._lowest_x = x
150 self._lowest_f = fx
152 return fx
154 def update_fun():
155 self.f = fun_wrapped(self.x)
157 self._update_fun_impl = update_fun
158 self._update_fun()
160 # Gradient evaluation
161 if callable(grad):
162 def grad_wrapped(x):
163 self.ngev += 1
164 return np.atleast_1d(grad(np.copy(x), *args))
166 def update_grad():
167 self.g = grad_wrapped(self.x)
169 elif grad in FD_METHODS:
170 def update_grad():
171 self._update_fun()
172 self.ngev += 1
173 self.g = approx_derivative(fun_wrapped, self.x, f0=self.f,
174 **finite_diff_options)
176 self._update_grad_impl = update_grad
177 self._update_grad()
179 # Hessian Evaluation
180 if callable(hess):
181 self.H = hess(np.copy(x0), *args)
182 self.H_updated = True
183 self.nhev += 1
185 if sps.issparse(self.H):
186 def hess_wrapped(x):
187 self.nhev += 1
188 return sps.csr_matrix(hess(np.copy(x), *args))
189 self.H = sps.csr_matrix(self.H)
191 elif isinstance(self.H, LinearOperator):
192 def hess_wrapped(x):
193 self.nhev += 1
194 return hess(np.copy(x), *args)
196 else:
197 def hess_wrapped(x):
198 self.nhev += 1
199 return np.atleast_2d(np.asarray(hess(np.copy(x), *args)))
200 self.H = np.atleast_2d(np.asarray(self.H))
202 def update_hess():
203 self.H = hess_wrapped(self.x)
205 elif hess in FD_METHODS:
206 def update_hess():
207 self._update_grad()
208 self.H = approx_derivative(grad_wrapped, self.x, f0=self.g,
209 **finite_diff_options)
210 return self.H
212 update_hess()
213 self.H_updated = True
214 elif isinstance(hess, HessianUpdateStrategy):
215 self.H = hess
216 self.H.initialize(self.n, 'hess')
217 self.H_updated = True
218 self.x_prev = None
219 self.g_prev = None
221 def update_hess():
222 self._update_grad()
223 self.H.update(self.x - self.x_prev, self.g - self.g_prev)
225 self._update_hess_impl = update_hess
227 if isinstance(hess, HessianUpdateStrategy):
228 def update_x(x):
229 self._update_grad()
230 self.x_prev = self.x
231 self.g_prev = self.g
232 # ensure that self.x is a copy of x. Don't store a reference
233 # otherwise the memoization doesn't work properly.
234 self.x = np.atleast_1d(x).astype(float)
235 self.f_updated = False
236 self.g_updated = False
237 self.H_updated = False
238 self._update_hess()
239 else:
240 def update_x(x):
241 # ensure that self.x is a copy of x. Don't store a reference
242 # otherwise the memoization doesn't work properly.
243 self.x = np.atleast_1d(x).astype(float)
244 self.f_updated = False
245 self.g_updated = False
246 self.H_updated = False
247 self._update_x_impl = update_x
249 def _update_fun(self):
250 if not self.f_updated:
251 self._update_fun_impl()
252 self.f_updated = True
254 def _update_grad(self):
255 if not self.g_updated:
256 self._update_grad_impl()
257 self.g_updated = True
259 def _update_hess(self):
260 if not self.H_updated:
261 self._update_hess_impl()
262 self.H_updated = True
264 def fun(self, x):
265 if not np.array_equal(x, self.x):
266 self._update_x_impl(x)
267 self._update_fun()
268 return self.f
270 def grad(self, x):
271 if not np.array_equal(x, self.x):
272 self._update_x_impl(x)
273 self._update_grad()
274 return self.g
276 def hess(self, x):
277 if not np.array_equal(x, self.x):
278 self._update_x_impl(x)
279 self._update_hess()
280 return self.H
282 def fun_and_grad(self, x):
283 if not np.array_equal(x, self.x):
284 self._update_x_impl(x)
285 self._update_fun()
286 self._update_grad()
287 return self.f, self.g
290class VectorFunction:
291 """Vector function and its derivatives.
293 This class defines a vector function F: R^n->R^m and methods for
294 computing or approximating its first and second derivatives.
296 Notes
297 -----
298 This class implements a memoization logic. There are methods `fun`,
299 `jac`, hess` and corresponding attributes `f`, `J` and `H`. The following
300 things should be considered:
302 1. Use only public methods `fun`, `jac` and `hess`.
303 2. After one of the methods is called, the corresponding attribute
304 will be set. However, a subsequent call with a different argument
305 of *any* of the methods may overwrite the attribute.
306 """
307 def __init__(self, fun, x0, jac, hess,
308 finite_diff_rel_step, finite_diff_jac_sparsity,
309 finite_diff_bounds, sparse_jacobian):
310 if not callable(jac) and jac not in FD_METHODS:
311 raise ValueError("`jac` must be either callable or one of {}."
312 .format(FD_METHODS))
314 if not (callable(hess) or hess in FD_METHODS
315 or isinstance(hess, HessianUpdateStrategy)):
316 raise ValueError("`hess` must be either callable,"
317 "HessianUpdateStrategy or one of {}."
318 .format(FD_METHODS))
320 if jac in FD_METHODS and hess in FD_METHODS:
321 raise ValueError("Whenever the Jacobian is estimated via "
322 "finite-differences, we require the Hessian to "
323 "be estimated using one of the quasi-Newton "
324 "strategies.")
326 self.x = np.atleast_1d(x0).astype(float)
327 self.n = self.x.size
328 self.nfev = 0
329 self.njev = 0
330 self.nhev = 0
331 self.f_updated = False
332 self.J_updated = False
333 self.H_updated = False
335 finite_diff_options = {}
336 if jac in FD_METHODS:
337 finite_diff_options["method"] = jac
338 finite_diff_options["rel_step"] = finite_diff_rel_step
339 if finite_diff_jac_sparsity is not None:
340 sparsity_groups = group_columns(finite_diff_jac_sparsity)
341 finite_diff_options["sparsity"] = (finite_diff_jac_sparsity,
342 sparsity_groups)
343 finite_diff_options["bounds"] = finite_diff_bounds
344 self.x_diff = np.copy(self.x)
345 if hess in FD_METHODS:
346 finite_diff_options["method"] = hess
347 finite_diff_options["rel_step"] = finite_diff_rel_step
348 finite_diff_options["as_linear_operator"] = True
349 self.x_diff = np.copy(self.x)
350 if jac in FD_METHODS and hess in FD_METHODS:
351 raise ValueError("Whenever the Jacobian is estimated via "
352 "finite-differences, we require the Hessian to "
353 "be estimated using one of the quasi-Newton "
354 "strategies.")
356 # Function evaluation
357 def fun_wrapped(x):
358 self.nfev += 1
359 return np.atleast_1d(fun(x))
361 def update_fun():
362 self.f = fun_wrapped(self.x)
364 self._update_fun_impl = update_fun
365 update_fun()
367 self.v = np.zeros_like(self.f)
368 self.m = self.v.size
370 # Jacobian Evaluation
371 if callable(jac):
372 self.J = jac(self.x)
373 self.J_updated = True
374 self.njev += 1
376 if (sparse_jacobian or
377 sparse_jacobian is None and sps.issparse(self.J)):
378 def jac_wrapped(x):
379 self.njev += 1
380 return sps.csr_matrix(jac(x))
381 self.J = sps.csr_matrix(self.J)
382 self.sparse_jacobian = True
384 elif sps.issparse(self.J):
385 def jac_wrapped(x):
386 self.njev += 1
387 return jac(x).toarray()
388 self.J = self.J.toarray()
389 self.sparse_jacobian = False
391 else:
392 def jac_wrapped(x):
393 self.njev += 1
394 return np.atleast_2d(jac(x))
395 self.J = np.atleast_2d(self.J)
396 self.sparse_jacobian = False
398 def update_jac():
399 self.J = jac_wrapped(self.x)
401 elif jac in FD_METHODS:
402 self.J = approx_derivative(fun_wrapped, self.x, f0=self.f,
403 **finite_diff_options)
404 self.J_updated = True
406 if (sparse_jacobian or
407 sparse_jacobian is None and sps.issparse(self.J)):
408 def update_jac():
409 self._update_fun()
410 self.J = sps.csr_matrix(
411 approx_derivative(fun_wrapped, self.x, f0=self.f,
412 **finite_diff_options))
413 self.J = sps.csr_matrix(self.J)
414 self.sparse_jacobian = True
416 elif sps.issparse(self.J):
417 def update_jac():
418 self._update_fun()
419 self.J = approx_derivative(fun_wrapped, self.x, f0=self.f,
420 **finite_diff_options).toarray()
421 self.J = self.J.toarray()
422 self.sparse_jacobian = False
424 else:
425 def update_jac():
426 self._update_fun()
427 self.J = np.atleast_2d(
428 approx_derivative(fun_wrapped, self.x, f0=self.f,
429 **finite_diff_options))
430 self.J = np.atleast_2d(self.J)
431 self.sparse_jacobian = False
433 self._update_jac_impl = update_jac
435 # Define Hessian
436 if callable(hess):
437 self.H = hess(self.x, self.v)
438 self.H_updated = True
439 self.nhev += 1
441 if sps.issparse(self.H):
442 def hess_wrapped(x, v):
443 self.nhev += 1
444 return sps.csr_matrix(hess(x, v))
445 self.H = sps.csr_matrix(self.H)
447 elif isinstance(self.H, LinearOperator):
448 def hess_wrapped(x, v):
449 self.nhev += 1
450 return hess(x, v)
452 else:
453 def hess_wrapped(x, v):
454 self.nhev += 1
455 return np.atleast_2d(np.asarray(hess(x, v)))
456 self.H = np.atleast_2d(np.asarray(self.H))
458 def update_hess():
459 self.H = hess_wrapped(self.x, self.v)
460 elif hess in FD_METHODS:
461 def jac_dot_v(x, v):
462 return jac_wrapped(x).T.dot(v)
464 def update_hess():
465 self._update_jac()
466 self.H = approx_derivative(jac_dot_v, self.x,
467 f0=self.J.T.dot(self.v),
468 args=(self.v,),
469 **finite_diff_options)
470 update_hess()
471 self.H_updated = True
472 elif isinstance(hess, HessianUpdateStrategy):
473 self.H = hess
474 self.H.initialize(self.n, 'hess')
475 self.H_updated = True
476 self.x_prev = None
477 self.J_prev = None
479 def update_hess():
480 self._update_jac()
481 # When v is updated before x was updated, then x_prev and
482 # J_prev are None and we need this check.
483 if self.x_prev is not None and self.J_prev is not None:
484 delta_x = self.x - self.x_prev
485 delta_g = self.J.T.dot(self.v) - self.J_prev.T.dot(self.v)
486 self.H.update(delta_x, delta_g)
488 self._update_hess_impl = update_hess
490 if isinstance(hess, HessianUpdateStrategy):
491 def update_x(x):
492 self._update_jac()
493 self.x_prev = self.x
494 self.J_prev = self.J
495 self.x = np.atleast_1d(x).astype(float)
496 self.f_updated = False
497 self.J_updated = False
498 self.H_updated = False
499 self._update_hess()
500 else:
501 def update_x(x):
502 self.x = np.atleast_1d(x).astype(float)
503 self.f_updated = False
504 self.J_updated = False
505 self.H_updated = False
507 self._update_x_impl = update_x
509 def _update_v(self, v):
510 if not np.array_equal(v, self.v):
511 self.v = v
512 self.H_updated = False
514 def _update_x(self, x):
515 if not np.array_equal(x, self.x):
516 self._update_x_impl(x)
518 def _update_fun(self):
519 if not self.f_updated:
520 self._update_fun_impl()
521 self.f_updated = True
523 def _update_jac(self):
524 if not self.J_updated:
525 self._update_jac_impl()
526 self.J_updated = True
528 def _update_hess(self):
529 if not self.H_updated:
530 self._update_hess_impl()
531 self.H_updated = True
533 def fun(self, x):
534 self._update_x(x)
535 self._update_fun()
536 return self.f
538 def jac(self, x):
539 self._update_x(x)
540 self._update_jac()
541 return self.J
543 def hess(self, x, v):
544 # v should be updated before x.
545 self._update_v(v)
546 self._update_x(x)
547 self._update_hess()
548 return self.H
551class LinearVectorFunction:
552 """Linear vector function and its derivatives.
554 Defines a linear function F = A x, where x is N-D vector and
555 A is m-by-n matrix. The Jacobian is constant and equals to A. The Hessian
556 is identically zero and it is returned as a csr matrix.
557 """
558 def __init__(self, A, x0, sparse_jacobian):
559 if sparse_jacobian or sparse_jacobian is None and sps.issparse(A):
560 self.J = sps.csr_matrix(A)
561 self.sparse_jacobian = True
562 elif sps.issparse(A):
563 self.J = A.toarray()
564 self.sparse_jacobian = False
565 else:
566 # np.asarray makes sure A is ndarray and not matrix
567 self.J = np.atleast_2d(np.asarray(A))
568 self.sparse_jacobian = False
570 self.m, self.n = self.J.shape
572 self.x = np.atleast_1d(x0).astype(float)
573 self.f = self.J.dot(self.x)
574 self.f_updated = True
576 self.v = np.zeros(self.m, dtype=float)
577 self.H = sps.csr_matrix((self.n, self.n))
579 def _update_x(self, x):
580 if not np.array_equal(x, self.x):
581 self.x = np.atleast_1d(x).astype(float)
582 self.f_updated = False
584 def fun(self, x):
585 self._update_x(x)
586 if not self.f_updated:
587 self.f = self.J.dot(x)
588 self.f_updated = True
589 return self.f
591 def jac(self, x):
592 self._update_x(x)
593 return self.J
595 def hess(self, x, v):
596 self._update_x(x)
597 self.v = v
598 return self.H
601class IdentityVectorFunction(LinearVectorFunction):
602 """Identity vector function and its derivatives.
604 The Jacobian is the identity matrix, returned as a dense array when
605 `sparse_jacobian=False` and as a csr matrix otherwise. The Hessian is
606 identically zero and it is returned as a csr matrix.
607 """
608 def __init__(self, x0, sparse_jacobian):
609 n = len(x0)
610 if sparse_jacobian or sparse_jacobian is None:
611 A = sps.eye(n, format='csr')
612 sparse_jacobian = True
613 else:
614 A = np.eye(n)
615 sparse_jacobian = False
616 super().__init__(A, x0, sparse_jacobian)