Coverage for /usr/lib/python3/dist-packages/scipy/optimize/_shgo.py: 9%
595 statements
« prev ^ index » next coverage.py v7.9.1, created at 2025-06-14 15:55 +0200
1"""shgo: The simplicial homology global optimisation algorithm."""
2from collections import namedtuple
3import time
4import logging
5import warnings
6import sys
8import numpy as np
10from scipy import spatial
11from scipy.optimize import OptimizeResult, minimize, Bounds
12from scipy.optimize._optimize import MemoizeJac
13from scipy.optimize._constraints import new_bounds_to_old
14from scipy.optimize._minimize import standardize_constraints
15from scipy._lib._util import _FunctionWrapper
17from scipy.optimize._shgo_lib._complex import Complex
19__all__ = ['shgo']
22def shgo(
23 func, bounds, args=(), constraints=None, n=100, iters=1, callback=None,
24 minimizer_kwargs=None, options=None, sampling_method='simplicial', *,
25 workers=1
26):
27 """
28 Finds the global minimum of a function using SHG optimization.
30 SHGO stands for "simplicial homology global optimization".
32 Parameters
33 ----------
34 func : callable
35 The objective function to be minimized. Must be in the form
36 ``f(x, *args)``, where ``x`` is the argument in the form of a 1-D array
37 and ``args`` is a tuple of any additional fixed parameters needed to
38 completely specify the function.
39 bounds : sequence or `Bounds`
40 Bounds for variables. There are two ways to specify the bounds:
42 1. Instance of `Bounds` class.
43 2. Sequence of ``(min, max)`` pairs for each element in `x`.
45 args : tuple, optional
46 Any additional fixed parameters needed to completely specify the
47 objective function.
48 constraints : {Constraint, dict} or List of {Constraint, dict}, optional
49 Constraints definition. Only for COBYLA, SLSQP and trust-constr.
50 See the tutorial [5]_ for further details on specifying constraints.
52 .. note::
54 Only COBYLA, SLSQP, and trust-constr local minimize methods
55 currently support constraint arguments. If the ``constraints``
56 sequence used in the local optimization problem is not defined in
57 ``minimizer_kwargs`` and a constrained method is used then the
58 global ``constraints`` will be used.
59 (Defining a ``constraints`` sequence in ``minimizer_kwargs``
60 means that ``constraints`` will not be added so if equality
61 constraints and so forth need to be added then the inequality
62 functions in ``constraints`` need to be added to
63 ``minimizer_kwargs`` too).
64 COBYLA only supports inequality constraints.
66 .. versionchanged:: 1.11.0
68 ``constraints`` accepts `NonlinearConstraint`, `LinearConstraint`.
70 n : int, optional
71 Number of sampling points used in the construction of the simplicial
72 complex. For the default ``simplicial`` sampling method 2**dim + 1
73 sampling points are generated instead of the default `n=100`. For all
74 other specified values `n` sampling points are generated. For
75 ``sobol``, ``halton`` and other arbitrary `sampling_methods` `n=100` or
76 another speciefied number of sampling points are generated.
77 iters : int, optional
78 Number of iterations used in the construction of the simplicial
79 complex. Default is 1.
80 callback : callable, optional
81 Called after each iteration, as ``callback(xk)``, where ``xk`` is the
82 current parameter vector.
83 minimizer_kwargs : dict, optional
84 Extra keyword arguments to be passed to the minimizer
85 ``scipy.optimize.minimize`` Some important options could be:
87 * method : str
88 The minimization method. If not given, chosen to be one of
89 BFGS, L-BFGS-B, SLSQP, depending on whether or not the
90 problem has constraints or bounds.
91 * args : tuple
92 Extra arguments passed to the objective function (``func``) and
93 its derivatives (Jacobian, Hessian).
94 * options : dict, optional
95 Note that by default the tolerance is specified as
96 ``{ftol: 1e-12}``
98 options : dict, optional
99 A dictionary of solver options. Many of the options specified for the
100 global routine are also passed to the scipy.optimize.minimize routine.
101 The options that are also passed to the local routine are marked with
102 "(L)".
104 Stopping criteria, the algorithm will terminate if any of the specified
105 criteria are met. However, the default algorithm does not require any
106 to be specified:
108 * maxfev : int (L)
109 Maximum number of function evaluations in the feasible domain.
110 (Note only methods that support this option will terminate
111 the routine at precisely exact specified value. Otherwise the
112 criterion will only terminate during a global iteration)
113 * f_min
114 Specify the minimum objective function value, if it is known.
115 * f_tol : float
116 Precision goal for the value of f in the stopping
117 criterion. Note that the global routine will also
118 terminate if a sampling point in the global routine is
119 within this tolerance.
120 * maxiter : int
121 Maximum number of iterations to perform.
122 * maxev : int
123 Maximum number of sampling evaluations to perform (includes
124 searching in infeasible points).
125 * maxtime : float
126 Maximum processing runtime allowed
127 * minhgrd : int
128 Minimum homology group rank differential. The homology group of the
129 objective function is calculated (approximately) during every
130 iteration. The rank of this group has a one-to-one correspondence
131 with the number of locally convex subdomains in the objective
132 function (after adequate sampling points each of these subdomains
133 contain a unique global minimum). If the difference in the hgr is 0
134 between iterations for ``maxhgrd`` specified iterations the
135 algorithm will terminate.
137 Objective function knowledge:
139 * symmetry : list or bool
140 Specify if the objective function contains symmetric variables.
141 The search space (and therefore performance) is decreased by up to
142 O(n!) times in the fully symmetric case. If `True` is specified
143 then all variables will be set symmetric to the first variable.
144 Default
145 is set to False.
147 E.g. f(x) = (x_1 + x_2 + x_3) + (x_4)**2 + (x_5)**2 + (x_6)**2
149 In this equation x_2 and x_3 are symmetric to x_1, while x_5 and
150 x_6 are symmetric to x_4, this can be specified to the solver as:
152 symmetry = [0, # Variable 1
153 0, # symmetric to variable 1
154 0, # symmetric to variable 1
155 3, # Variable 4
156 3, # symmetric to variable 4
157 3, # symmetric to variable 4
158 ]
160 * jac : bool or callable, optional
161 Jacobian (gradient) of objective function. Only for CG, BFGS,
162 Newton-CG, L-BFGS-B, TNC, SLSQP, dogleg, trust-ncg. If ``jac`` is a
163 boolean and is True, ``fun`` is assumed to return the gradient
164 along with the objective function. If False, the gradient will be
165 estimated numerically. ``jac`` can also be a callable returning the
166 gradient of the objective. In this case, it must accept the same
167 arguments as ``fun``. (Passed to `scipy.optimize.minmize`
168 automatically)
170 * hess, hessp : callable, optional
171 Hessian (matrix of second-order derivatives) of objective function
172 or Hessian of objective function times an arbitrary vector p.
173 Only for Newton-CG, dogleg, trust-ncg. Only one of ``hessp`` or
174 ``hess`` needs to be given. If ``hess`` is provided, then
175 ``hessp`` will be ignored. If neither ``hess`` nor ``hessp`` is
176 provided, then the Hessian product will be approximated using
177 finite differences on ``jac``. ``hessp`` must compute the Hessian
178 times an arbitrary vector. (Passed to `scipy.optimize.minmize`
179 automatically)
181 Algorithm settings:
183 * minimize_every_iter : bool
184 If True then promising global sampling points will be passed to a
185 local minimization routine every iteration. If True then only the
186 final minimizer pool will be run. Defaults to True.
187 * local_iter : int
188 Only evaluate a few of the best minimizer pool candidates every
189 iteration. If False all potential points are passed to the local
190 minimization routine.
191 * infty_constraints : bool
192 If True then any sampling points generated which are outside will
193 the feasible domain will be saved and given an objective function
194 value of ``inf``. If False then these points will be discarded.
195 Using this functionality could lead to higher performance with
196 respect to function evaluations before the global minimum is found,
197 specifying False will use less memory at the cost of a slight
198 decrease in performance. Defaults to True.
200 Feedback:
202 * disp : bool (L)
203 Set to True to print convergence messages.
205 sampling_method : str or function, optional
206 Current built in sampling method options are ``halton``, ``sobol`` and
207 ``simplicial``. The default ``simplicial`` provides
208 the theoretical guarantee of convergence to the global minimum in
209 finite time. ``halton`` and ``sobol`` method are faster in terms of
210 sampling point generation at the cost of the loss of
211 guaranteed convergence. It is more appropriate for most "easier"
212 problems where the convergence is relatively fast.
213 User defined sampling functions must accept two arguments of ``n``
214 sampling points of dimension ``dim`` per call and output an array of
215 sampling points with shape `n x dim`.
217 workers : int or map-like callable, optional
218 Sample and run the local serial minimizations in parallel.
219 Supply -1 to use all available CPU cores, or an int to use
220 that many Processes (uses `multiprocessing.Pool <multiprocessing>`).
222 Alternatively supply a map-like callable, such as
223 `multiprocessing.Pool.map` for parallel evaluation.
224 This evaluation is carried out as ``workers(func, iterable)``.
225 Requires that `func` be pickleable.
227 .. versionadded:: 1.11.0
229 Returns
230 -------
231 res : OptimizeResult
232 The optimization result represented as a `OptimizeResult` object.
233 Important attributes are:
234 ``x`` the solution array corresponding to the global minimum,
235 ``fun`` the function output at the global solution,
236 ``xl`` an ordered list of local minima solutions,
237 ``funl`` the function output at the corresponding local solutions,
238 ``success`` a Boolean flag indicating if the optimizer exited
239 successfully,
240 ``message`` which describes the cause of the termination,
241 ``nfev`` the total number of objective function evaluations including
242 the sampling calls,
243 ``nlfev`` the total number of objective function evaluations
244 culminating from all local search optimizations,
245 ``nit`` number of iterations performed by the global routine.
247 Notes
248 -----
249 Global optimization using simplicial homology global optimization [1]_.
250 Appropriate for solving general purpose NLP and blackbox optimization
251 problems to global optimality (low-dimensional problems).
253 In general, the optimization problems are of the form::
255 minimize f(x) subject to
257 g_i(x) >= 0, i = 1,...,m
258 h_j(x) = 0, j = 1,...,p
260 where x is a vector of one or more variables. ``f(x)`` is the objective
261 function ``R^n -> R``, ``g_i(x)`` are the inequality constraints, and
262 ``h_j(x)`` are the equality constraints.
264 Optionally, the lower and upper bounds for each element in x can also be
265 specified using the `bounds` argument.
267 While most of the theoretical advantages of SHGO are only proven for when
268 ``f(x)`` is a Lipschitz smooth function, the algorithm is also proven to
269 converge to the global optimum for the more general case where ``f(x)`` is
270 non-continuous, non-convex and non-smooth, if the default sampling method
271 is used [1]_.
273 The local search method may be specified using the ``minimizer_kwargs``
274 parameter which is passed on to ``scipy.optimize.minimize``. By default,
275 the ``SLSQP`` method is used. In general, it is recommended to use the
276 ``SLSQP`` or ``COBYLA`` local minimization if inequality constraints
277 are defined for the problem since the other methods do not use constraints.
279 The ``halton`` and ``sobol`` method points are generated using
280 `scipy.stats.qmc`. Any other QMC method could be used.
282 References
283 ----------
284 .. [1] Endres, SC, Sandrock, C, Focke, WW (2018) "A simplicial homology
285 algorithm for lipschitz optimisation", Journal of Global
286 Optimization.
287 .. [2] Joe, SW and Kuo, FY (2008) "Constructing Sobol' sequences with
288 better two-dimensional projections", SIAM J. Sci. Comput. 30,
289 2635-2654.
290 .. [3] Hock, W and Schittkowski, K (1981) "Test examples for nonlinear
291 programming codes", Lecture Notes in Economics and Mathematical
292 Systems, 187. Springer-Verlag, New York.
293 http://www.ai7.uni-bayreuth.de/test_problem_coll.pdf
294 .. [4] Wales, DJ (2015) "Perspective: Insight into reaction coordinates and
295 dynamics from the potential energy landscape",
296 Journal of Chemical Physics, 142(13), 2015.
297 .. [5] https://docs.scipy.org/doc/scipy/tutorial/optimize.html#constrained-minimization-of-multivariate-scalar-functions-minimize
299 Examples
300 --------
301 First consider the problem of minimizing the Rosenbrock function, `rosen`:
303 >>> from scipy.optimize import rosen, shgo
304 >>> bounds = [(0,2), (0, 2), (0, 2), (0, 2), (0, 2)]
305 >>> result = shgo(rosen, bounds)
306 >>> result.x, result.fun
307 (array([1., 1., 1., 1., 1.]), 2.920392374190081e-18)
309 Note that bounds determine the dimensionality of the objective
310 function and is therefore a required input, however you can specify
311 empty bounds using ``None`` or objects like ``np.inf`` which will be
312 converted to large float numbers.
314 >>> bounds = [(None, None), ]*4
315 >>> result = shgo(rosen, bounds)
316 >>> result.x
317 array([0.99999851, 0.99999704, 0.99999411, 0.9999882 ])
319 Next, we consider the Eggholder function, a problem with several local
320 minima and one global minimum. We will demonstrate the use of arguments and
321 the capabilities of `shgo`.
322 (https://en.wikipedia.org/wiki/Test_functions_for_optimization)
324 >>> import numpy as np
325 >>> def eggholder(x):
326 ... return (-(x[1] + 47.0)
327 ... * np.sin(np.sqrt(abs(x[0]/2.0 + (x[1] + 47.0))))
328 ... - x[0] * np.sin(np.sqrt(abs(x[0] - (x[1] + 47.0))))
329 ... )
330 ...
331 >>> bounds = [(-512, 512), (-512, 512)]
333 `shgo` has built-in low discrepancy sampling sequences. First, we will
334 input 64 initial sampling points of the *Sobol'* sequence:
336 >>> result = shgo(eggholder, bounds, n=64, sampling_method='sobol')
337 >>> result.x, result.fun
338 (array([512. , 404.23180824]), -959.6406627208397)
340 `shgo` also has a return for any other local minima that was found, these
341 can be called using:
343 >>> result.xl
344 array([[ 512. , 404.23180824],
345 [ 283.0759062 , -487.12565635],
346 [-294.66820039, -462.01964031],
347 [-105.87688911, 423.15323845],
348 [-242.97926 , 274.38030925],
349 [-506.25823477, 6.3131022 ],
350 [-408.71980731, -156.10116949],
351 [ 150.23207937, 301.31376595],
352 [ 91.00920901, -391.283763 ],
353 [ 202.89662724, -269.38043241],
354 [ 361.66623976, -106.96493868],
355 [-219.40612786, -244.06020508]])
357 >>> result.funl
358 array([-959.64066272, -718.16745962, -704.80659592, -565.99778097,
359 -559.78685655, -557.36868733, -507.87385942, -493.9605115 ,
360 -426.48799655, -421.15571437, -419.31194957, -410.98477763])
362 These results are useful in applications where there are many global minima
363 and the values of other global minima are desired or where the local minima
364 can provide insight into the system (for example morphologies
365 in physical chemistry [4]_).
367 If we want to find a larger number of local minima, we can increase the
368 number of sampling points or the number of iterations. We'll increase the
369 number of sampling points to 64 and the number of iterations from the
370 default of 1 to 3. Using ``simplicial`` this would have given us
371 64 x 3 = 192 initial sampling points.
373 >>> result_2 = shgo(eggholder,
374 ... bounds, n=64, iters=3, sampling_method='sobol')
375 >>> len(result.xl), len(result_2.xl)
376 (12, 23)
378 Note the difference between, e.g., ``n=192, iters=1`` and ``n=64,
379 iters=3``.
380 In the first case the promising points contained in the minimiser pool
381 are processed only once. In the latter case it is processed every 64
382 sampling points for a total of 3 times.
384 To demonstrate solving problems with non-linear constraints consider the
385 following example from Hock and Schittkowski problem 73 (cattle-feed)
386 [3]_::
388 minimize: f = 24.55 * x_1 + 26.75 * x_2 + 39 * x_3 + 40.50 * x_4
390 subject to: 2.3 * x_1 + 5.6 * x_2 + 11.1 * x_3 + 1.3 * x_4 - 5 >= 0,
392 12 * x_1 + 11.9 * x_2 + 41.8 * x_3 + 52.1 * x_4 - 21
393 -1.645 * sqrt(0.28 * x_1**2 + 0.19 * x_2**2 +
394 20.5 * x_3**2 + 0.62 * x_4**2) >= 0,
396 x_1 + x_2 + x_3 + x_4 - 1 == 0,
398 1 >= x_i >= 0 for all i
400 The approximate answer given in [3]_ is::
402 f([0.6355216, -0.12e-11, 0.3127019, 0.05177655]) = 29.894378
404 >>> def f(x): # (cattle-feed)
405 ... return 24.55*x[0] + 26.75*x[1] + 39*x[2] + 40.50*x[3]
406 ...
407 >>> def g1(x):
408 ... return 2.3*x[0] + 5.6*x[1] + 11.1*x[2] + 1.3*x[3] - 5 # >=0
409 ...
410 >>> def g2(x):
411 ... return (12*x[0] + 11.9*x[1] +41.8*x[2] + 52.1*x[3] - 21
412 ... - 1.645 * np.sqrt(0.28*x[0]**2 + 0.19*x[1]**2
413 ... + 20.5*x[2]**2 + 0.62*x[3]**2)
414 ... ) # >=0
415 ...
416 >>> def h1(x):
417 ... return x[0] + x[1] + x[2] + x[3] - 1 # == 0
418 ...
419 >>> cons = ({'type': 'ineq', 'fun': g1},
420 ... {'type': 'ineq', 'fun': g2},
421 ... {'type': 'eq', 'fun': h1})
422 >>> bounds = [(0, 1.0),]*4
423 >>> res = shgo(f, bounds, n=150, constraints=cons)
424 >>> res
425 message: Optimization terminated successfully.
426 success: True
427 fun: 29.894378159142136
428 funl: [ 2.989e+01]
429 x: [ 6.355e-01 1.137e-13 3.127e-01 5.178e-02]
430 xl: [[ 6.355e-01 1.137e-13 3.127e-01 5.178e-02]]
431 nit: 1
432 nfev: 142
433 nlfev: 35
434 nljev: 5
435 nlhev: 0
437 >>> g1(res.x), g2(res.x), h1(res.x)
438 (-5.062616992290714e-14, -2.9594104944408173e-12, 0.0)
440 """
441 # if necessary, convert bounds class to old bounds
442 if isinstance(bounds, Bounds):
443 bounds = new_bounds_to_old(bounds.lb, bounds.ub, len(bounds.lb))
445 # Initiate SHGO class
446 # use in context manager to make sure that any parallelization
447 # resources are freed.
448 with SHGO(func, bounds, args=args, constraints=constraints, n=n,
449 iters=iters, callback=callback,
450 minimizer_kwargs=minimizer_kwargs,
451 options=options, sampling_method=sampling_method,
452 workers=workers) as shc:
453 # Run the algorithm, process results and test success
454 shc.iterate_all()
456 if not shc.break_routine:
457 if shc.disp:
458 logging.info("Successfully completed construction of complex.")
460 # Test post iterations success
461 if len(shc.LMC.xl_maps) == 0:
462 # If sampling failed to find pool, return lowest sampled point
463 # with a warning
464 shc.find_lowest_vertex()
465 shc.break_routine = True
466 shc.fail_routine(mes="Failed to find a feasible minimizer point. "
467 "Lowest sampling point = {}".format(shc.f_lowest))
468 shc.res.fun = shc.f_lowest
469 shc.res.x = shc.x_lowest
470 shc.res.nfev = shc.fn
471 shc.res.tnev = shc.n_sampled
472 else:
473 # Test that the optimal solutions do not violate any constraints
474 pass # TODO
476 # Confirm the routine ran successfully
477 if not shc.break_routine:
478 shc.res.message = 'Optimization terminated successfully.'
479 shc.res.success = True
481 # Return the final results
482 return shc.res
485class SHGO:
486 def __init__(self, func, bounds, args=(), constraints=None, n=None,
487 iters=None, callback=None, minimizer_kwargs=None,
488 options=None, sampling_method='simplicial', workers=1):
489 from scipy.stats import qmc
490 # Input checks
491 methods = ['halton', 'sobol', 'simplicial']
492 if isinstance(sampling_method, str) and sampling_method not in methods:
493 raise ValueError(("Unknown sampling_method specified."
494 " Valid methods: {}").format(', '.join(methods)))
496 # Split obj func if given with Jac
497 try:
498 if ((minimizer_kwargs['jac'] is True) and
499 (not callable(minimizer_kwargs['jac']))):
500 self.func = MemoizeJac(func)
501 jac = self.func.derivative
502 minimizer_kwargs['jac'] = jac
503 func = self.func # .fun
504 else:
505 self.func = func # Normal definition of objective function
506 except (TypeError, KeyError):
507 self.func = func # Normal definition of objective function
509 # Initiate class
510 self.func = _FunctionWrapper(func, args)
511 self.bounds = bounds
512 self.args = args
513 self.callback = callback
515 # Bounds
516 abound = np.array(bounds, float)
517 self.dim = np.shape(abound)[0] # Dimensionality of problem
519 # Set none finite values to large floats
520 infind = ~np.isfinite(abound)
521 abound[infind[:, 0], 0] = -1e50
522 abound[infind[:, 1], 1] = 1e50
524 # Check if bounds are correctly specified
525 bnderr = abound[:, 0] > abound[:, 1]
526 if bnderr.any():
527 raise ValueError('Error: lb > ub in bounds {}.'
528 .format(', '.join(str(b) for b in bnderr)))
530 self.bounds = abound
532 # Constraints
533 # Process constraint dict sequence:
534 self.constraints = constraints
535 if constraints is not None:
536 self.min_cons = constraints
537 self.g_cons = []
538 self.g_args = []
540 # shgo internals deals with old-style constraints
541 # self.constraints is used to create Complex, so need
542 # to be stored internally in old-style.
543 # `minimize` takes care of normalising these constraints
544 # for slsqp/cobyla/trust-constr.
545 self.constraints = standardize_constraints(
546 constraints,
547 np.empty(self.dim, float),
548 'old'
549 )
550 for cons in self.constraints:
551 if cons['type'] in ('ineq'):
552 self.g_cons.append(cons['fun'])
553 try:
554 self.g_args.append(cons['args'])
555 except KeyError:
556 self.g_args.append(())
557 self.g_cons = tuple(self.g_cons)
558 self.g_args = tuple(self.g_args)
559 else:
560 self.g_cons = None
561 self.g_args = None
563 # Define local minimization keyword arguments
564 # Start with defaults
565 self.minimizer_kwargs = {'method': 'SLSQP',
566 'bounds': self.bounds,
567 'options': {},
568 'callback': self.callback
569 }
570 if minimizer_kwargs is not None:
571 # Overwrite with supplied values
572 self.minimizer_kwargs.update(minimizer_kwargs)
574 else:
575 self.minimizer_kwargs['options'] = {'ftol': 1e-12}
577 if (
578 self.minimizer_kwargs['method'].lower() in ('slsqp', 'cobyla', 'trust-constr') and
579 (
580 minimizer_kwargs is not None and
581 'constraints' not in minimizer_kwargs and
582 constraints is not None
583 ) or
584 (self.g_cons is not None)
585 ):
586 self.minimizer_kwargs['constraints'] = self.min_cons
588 # Process options dict
589 if options is not None:
590 self.init_options(options)
591 else: # Default settings:
592 self.f_min_true = None
593 self.minimize_every_iter = True
595 # Algorithm limits
596 self.maxiter = None
597 self.maxfev = None
598 self.maxev = None
599 self.maxtime = None
600 self.f_min_true = None
601 self.minhgrd = None
603 # Objective function knowledge
604 self.symmetry = None
606 # Algorithm functionality
607 self.infty_cons_sampl = True
608 self.local_iter = False
610 # Feedback
611 self.disp = False
613 # Remove unknown arguments in self.minimizer_kwargs
614 # Start with arguments all the solvers have in common
615 self.min_solver_args = ['fun', 'x0', 'args',
616 'callback', 'options', 'method']
617 # then add the ones unique to specific solvers
618 solver_args = {
619 '_custom': ['jac', 'hess', 'hessp', 'bounds', 'constraints'],
620 'nelder-mead': [],
621 'powell': [],
622 'cg': ['jac'],
623 'bfgs': ['jac'],
624 'newton-cg': ['jac', 'hess', 'hessp'],
625 'l-bfgs-b': ['jac', 'bounds'],
626 'tnc': ['jac', 'bounds'],
627 'cobyla': ['constraints', 'catol'],
628 'slsqp': ['jac', 'bounds', 'constraints'],
629 'dogleg': ['jac', 'hess'],
630 'trust-ncg': ['jac', 'hess', 'hessp'],
631 'trust-krylov': ['jac', 'hess', 'hessp'],
632 'trust-exact': ['jac', 'hess'],
633 'trust-constr': ['jac', 'hess', 'hessp', 'constraints'],
634 }
635 method = self.minimizer_kwargs['method']
636 self.min_solver_args += solver_args[method.lower()]
638 # Only retain the known arguments
639 def _restrict_to_keys(dictionary, goodkeys):
640 """Remove keys from dictionary if not in goodkeys - inplace"""
641 existingkeys = set(dictionary)
642 for key in existingkeys - set(goodkeys):
643 dictionary.pop(key, None)
645 _restrict_to_keys(self.minimizer_kwargs, self.min_solver_args)
646 _restrict_to_keys(self.minimizer_kwargs['options'],
647 self.min_solver_args + ['ftol'])
649 # Algorithm controls
650 # Global controls
651 self.stop_global = False # Used in the stopping_criteria method
652 self.break_routine = False # Break the algorithm globally
653 self.iters = iters # Iterations to be ran
654 self.iters_done = 0 # Iterations completed
655 self.n = n # Sampling points per iteration
656 self.nc = 0 # n # Sampling points to sample in current iteration
657 self.n_prc = 0 # Processed points (used to track Delaunay iters)
658 self.n_sampled = 0 # To track no. of sampling points already generated
659 self.fn = 0 # Number of feasible sampling points evaluations performed
660 self.hgr = 0 # Homology group rank
661 # Initially attempt to build the triangulation incrementally:
662 self.qhull_incremental = True
664 # Default settings if no sampling criteria.
665 if (self.n is None) and (self.iters is None) \
666 and (sampling_method == 'simplicial'):
667 self.n = 2 ** self.dim + 1
668 self.nc = 0 # self.n
669 if self.iters is None:
670 self.iters = 1
671 if (self.n is None) and not (sampling_method == 'simplicial'):
672 self.n = self.n = 100
673 self.nc = 0 # self.n
674 if (self.n == 100) and (sampling_method == 'simplicial'):
675 self.n = 2 ** self.dim + 1
677 if not ((self.maxiter is None) and (self.maxfev is None) and (
678 self.maxev is None)
679 and (self.minhgrd is None) and (self.f_min_true is None)):
680 self.iters = None
682 # Set complex construction mode based on a provided stopping criteria:
683 # Initialise sampling Complex and function cache
684 # Note that sfield_args=() since args are already wrapped in self.func
685 # using the_FunctionWrapper class.
686 self.HC = Complex(dim=self.dim, domain=self.bounds,
687 sfield=self.func, sfield_args=(),
688 symmetry=self.symmetry,
689 constraints=self.constraints,
690 workers=workers)
692 # Choose complex constructor
693 if sampling_method == 'simplicial':
694 self.iterate_complex = self.iterate_hypercube
695 self.sampling_method = sampling_method
697 elif sampling_method in ['halton', 'sobol'] or \
698 not isinstance(sampling_method, str):
699 self.iterate_complex = self.iterate_delaunay
700 # Sampling method used
701 if sampling_method in ['halton', 'sobol']:
702 if sampling_method == 'sobol':
703 self.n = int(2 ** np.ceil(np.log2(self.n)))
704 # self.n #TODO: Should always be self.n, this is
705 # unacceptable for shgo, check that nfev behaves as
706 # expected.
707 self.nc = 0
708 self.sampling_method = 'sobol'
709 self.qmc_engine = qmc.Sobol(d=self.dim, scramble=False,
710 seed=0)
711 else:
712 self.sampling_method = 'halton'
713 self.qmc_engine = qmc.Halton(d=self.dim, scramble=True,
714 seed=0)
716 def sampling_method(n, d):
717 return self.qmc_engine.random(n)
719 else:
720 # A user defined sampling method:
721 self.sampling_method = 'custom'
723 self.sampling = self.sampling_custom
724 self.sampling_function = sampling_method # F(n, d)
726 # Local controls
727 self.stop_l_iter = False # Local minimisation iterations
728 self.stop_complex_iter = False # Sampling iterations
730 # Initiate storage objects used in algorithm classes
731 self.minimizer_pool = []
733 # Cache of local minimizers mapped
734 self.LMC = LMapCache()
736 # Initialize return object
737 self.res = OptimizeResult() # scipy.optimize.OptimizeResult object
738 self.res.nfev = 0 # Includes each sampling point as func evaluation
739 self.res.nlfev = 0 # Local function evals for all minimisers
740 self.res.nljev = 0 # Local Jacobian evals for all minimisers
741 self.res.nlhev = 0 # Local Hessian evals for all minimisers
743 # Initiation aids
744 def init_options(self, options):
745 """
746 Initiates the options.
748 Can also be useful to change parameters after class initiation.
750 Parameters
751 ----------
752 options : dict
754 Returns
755 -------
756 None
758 """
759 # Update 'options' dict passed to optimize.minimize
760 # Do this first so we don't mutate `options` below.
761 self.minimizer_kwargs['options'].update(options)
763 # Ensure that 'jac', 'hess', and 'hessp' are passed directly to
764 # `minimize` as keywords, not as part of its 'options' dictionary.
765 for opt in ['jac', 'hess', 'hessp']:
766 if opt in self.minimizer_kwargs['options']:
767 self.minimizer_kwargs[opt] = (
768 self.minimizer_kwargs['options'].pop(opt))
770 # Default settings:
771 self.minimize_every_iter = options.get('minimize_every_iter', True)
773 # Algorithm limits
774 # Maximum number of iterations to perform.
775 self.maxiter = options.get('maxiter', None)
776 # Maximum number of function evaluations in the feasible domain
777 self.maxfev = options.get('maxfev', None)
778 # Maximum number of sampling evaluations (includes searching in
779 # infeasible points
780 self.maxev = options.get('maxev', None)
781 # Maximum processing runtime allowed
782 self.init = time.time()
783 self.maxtime = options.get('maxtime', None)
784 if 'f_min' in options:
785 # Specify the minimum objective function value, if it is known.
786 self.f_min_true = options['f_min']
787 self.f_tol = options.get('f_tol', 1e-4)
788 else:
789 self.f_min_true = None
791 self.minhgrd = options.get('minhgrd', None)
793 # Objective function knowledge
794 self.symmetry = options.get('symmetry', False)
795 if self.symmetry:
796 self.symmetry = [0, ]*len(self.bounds)
797 else:
798 self.symmetry = None
799 # Algorithm functionality
800 # Only evaluate a few of the best candiates
801 self.local_iter = options.get('local_iter', False)
802 self.infty_cons_sampl = options.get('infty_constraints', True)
804 # Feedback
805 self.disp = options.get('disp', False)
807 def __enter__(self):
808 return self
810 def __exit__(self, *args):
811 return self.HC.V._mapwrapper.__exit__(*args)
813 # Iteration properties
814 # Main construction loop:
815 def iterate_all(self):
816 """
817 Construct for `iters` iterations.
819 If uniform sampling is used, every iteration adds 'n' sampling points.
821 Iterations if a stopping criteria (e.g., sampling points or
822 processing time) has been met.
824 """
825 if self.disp:
826 logging.info('Splitting first generation')
828 while not self.stop_global:
829 if self.break_routine:
830 break
831 # Iterate complex, process minimisers
832 self.iterate()
833 self.stopping_criteria()
835 # Build minimiser pool
836 # Final iteration only needed if pools weren't minimised every
837 # iteration
838 if not self.minimize_every_iter:
839 if not self.break_routine:
840 self.find_minima()
842 self.res.nit = self.iters_done # + 1
843 self.fn = self.HC.V.nfev
845 def find_minima(self):
846 """
847 Construct the minimizer pool, map the minimizers to local minima
848 and sort the results into a global return object.
849 """
850 if self.disp:
851 logging.info('Searching for minimizer pool...')
853 self.minimizers()
855 if len(self.X_min) != 0:
856 # Minimize the pool of minimizers with local minimization methods
857 # Note that if Options['local_iter'] is an `int` instead of default
858 # value False then only that number of candidates will be minimized
859 self.minimise_pool(self.local_iter)
860 # Sort results and build the global return object
861 self.sort_result()
863 # Lowest values used to report in case of failures
864 self.f_lowest = self.res.fun
865 self.x_lowest = self.res.x
866 else:
867 self.find_lowest_vertex()
869 if self.disp:
870 logging.info(f"Minimiser pool = SHGO.X_min = {self.X_min}")
872 def find_lowest_vertex(self):
873 # Find the lowest objective function value on one of
874 # the vertices of the simplicial complex
875 self.f_lowest = np.inf
876 for x in self.HC.V.cache:
877 if self.HC.V[x].f < self.f_lowest:
878 if self.disp:
879 logging.info(f'self.HC.V[x].f = {self.HC.V[x].f}')
880 self.f_lowest = self.HC.V[x].f
881 self.x_lowest = self.HC.V[x].x_a
882 for lmc in self.LMC.cache:
883 if self.LMC[lmc].f_min < self.f_lowest:
884 self.f_lowest = self.LMC[lmc].f_min
885 self.x_lowest = self.LMC[lmc].x_l
887 if self.f_lowest == np.inf: # no feasible point
888 self.f_lowest = None
889 self.x_lowest = None
891 # Stopping criteria functions:
892 def finite_iterations(self):
893 mi = min(x for x in [self.iters, self.maxiter] if x is not None)
894 if self.disp:
895 logging.info(f'Iterations done = {self.iters_done} / {mi}')
896 if self.iters is not None:
897 if self.iters_done >= (self.iters):
898 self.stop_global = True
900 if self.maxiter is not None: # Stop for infeasible sampling
901 if self.iters_done >= (self.maxiter):
902 self.stop_global = True
903 return self.stop_global
905 def finite_fev(self):
906 # Finite function evals in the feasible domain
907 if self.disp:
908 logging.info(f'Function evaluations done = {self.fn} / {self.maxfev}')
909 if self.fn >= self.maxfev:
910 self.stop_global = True
911 return self.stop_global
913 def finite_ev(self):
914 # Finite evaluations including infeasible sampling points
915 if self.disp:
916 logging.info(f'Sampling evaluations done = {self.n_sampled} '
917 f'/ {self.maxev}')
918 if self.n_sampled >= self.maxev:
919 self.stop_global = True
921 def finite_time(self):
922 if self.disp:
923 logging.info(f'Time elapsed = {time.time() - self.init} '
924 f'/ {self.maxtime}')
925 if (time.time() - self.init) >= self.maxtime:
926 self.stop_global = True
928 def finite_precision(self):
929 """
930 Stop the algorithm if the final function value is known
932 Specify in options (with ``self.f_min_true = options['f_min']``)
933 and the tolerance with ``f_tol = options['f_tol']``
934 """
935 # If no minimizer has been found use the lowest sampling value
936 self.find_lowest_vertex()
937 if self.disp:
938 logging.info(f'Lowest function evaluation = {self.f_lowest}')
939 logging.info(f'Specified minimum = {self.f_min_true}')
940 # If no feasible point was return from test
941 if self.f_lowest is None:
942 return self.stop_global
944 # Function to stop algorithm at specified percentage error:
945 if self.f_min_true == 0.0:
946 if self.f_lowest <= self.f_tol:
947 self.stop_global = True
948 else:
949 pe = (self.f_lowest - self.f_min_true) / abs(self.f_min_true)
950 if self.f_lowest <= self.f_min_true:
951 self.stop_global = True
952 # 2if (pe - self.f_tol) <= abs(1.0 / abs(self.f_min_true)):
953 if abs(pe) >= 2 * self.f_tol:
954 warnings.warn("A much lower value than expected f* =" +
955 f" {self.f_min_true} than" +
956 " the was found f_lowest =" +
957 f"{self.f_lowest} ")
958 if pe <= self.f_tol:
959 self.stop_global = True
961 return self.stop_global
963 def finite_homology_growth(self):
964 """
965 Stop the algorithm if homology group rank did not grow in iteration.
966 """
967 if self.LMC.size == 0:
968 return # pass on no reason to stop yet.
969 self.hgrd = self.LMC.size - self.hgr
971 self.hgr = self.LMC.size
972 if self.hgrd <= self.minhgrd:
973 self.stop_global = True
974 if self.disp:
975 logging.info(f'Current homology growth = {self.hgrd} '
976 f' (minimum growth = {self.minhgrd})')
977 return self.stop_global
979 def stopping_criteria(self):
980 """
981 Various stopping criteria ran every iteration
983 Returns
984 -------
985 stop : bool
986 """
987 if self.maxiter is not None:
988 self.finite_iterations()
989 if self.iters is not None:
990 self.finite_iterations()
991 if self.maxfev is not None:
992 self.finite_fev()
993 if self.maxev is not None:
994 self.finite_ev()
995 if self.maxtime is not None:
996 self.finite_time()
997 if self.f_min_true is not None:
998 self.finite_precision()
999 if self.minhgrd is not None:
1000 self.finite_homology_growth()
1001 return self.stop_global
1003 def iterate(self):
1004 self.iterate_complex()
1006 # Build minimizer pool
1007 if self.minimize_every_iter:
1008 if not self.break_routine:
1009 self.find_minima() # Process minimizer pool
1011 # Algorithm updates
1012 self.iters_done += 1
1014 def iterate_hypercube(self):
1015 """
1016 Iterate a subdivision of the complex
1018 Note: called with ``self.iterate_complex()`` after class initiation
1019 """
1020 # Iterate the complex
1021 if self.disp:
1022 logging.info('Constructing and refining simplicial complex graph '
1023 'structure')
1024 if self.n is None:
1025 self.HC.refine_all()
1026 self.n_sampled = self.HC.V.size() # nevs counted
1027 else:
1028 self.HC.refine(self.n)
1029 self.n_sampled += self.n
1031 if self.disp:
1032 logging.info('Triangulation completed, evaluating all contraints '
1033 'and objective function values.')
1035 # Readd minimisers to complex
1036 if len(self.LMC.xl_maps) > 0:
1037 for xl in self.LMC.cache:
1038 v = self.HC.V[xl]
1039 v_near = v.star()
1040 for v in v.nn:
1041 v_near = v_near.union(v.nn)
1042 # Reconnect vertices to complex
1043 # if self.HC.connect_vertex_non_symm(tuple(self.LMC[xl].x_l),
1044 # near=v_near):
1045 # continue
1046 # else:
1047 # If failure to find in v_near, then search all vertices
1048 # (very expensive operation:
1049 # self.HC.connect_vertex_non_symm(tuple(self.LMC[xl].x_l)
1050 # )
1052 # Evaluate all constraints and functions
1053 self.HC.V.process_pools()
1054 if self.disp:
1055 logging.info('Evaluations completed.')
1057 # feasible sampling points counted by the triangulation.py routines
1058 self.fn = self.HC.V.nfev
1059 return
1061 def iterate_delaunay(self):
1062 """
1063 Build a complex of Delaunay triangulated points
1065 Note: called with ``self.iterate_complex()`` after class initiation
1066 """
1067 self.nc += self.n
1068 self.sampled_surface(infty_cons_sampl=self.infty_cons_sampl)
1070 # Add sampled points to a triangulation, construct self.Tri
1071 if self.disp:
1072 logging.info(f'self.n = {self.n}')
1073 logging.info(f'self.nc = {self.nc}')
1074 logging.info('Constructing and refining simplicial complex graph '
1075 'structure from sampling points.')
1077 if self.dim < 2:
1078 self.Ind_sorted = np.argsort(self.C, axis=0)
1079 self.Ind_sorted = self.Ind_sorted.flatten()
1080 tris = []
1081 for ind, ind_s in enumerate(self.Ind_sorted):
1082 if ind > 0:
1083 tris.append(self.Ind_sorted[ind - 1:ind + 1])
1085 tris = np.array(tris)
1086 # Store 1D triangulation:
1087 self.Tri = namedtuple('Tri', ['points', 'simplices'])(self.C, tris)
1088 self.points = {}
1089 else:
1090 if self.C.shape[0] > self.dim + 1: # Ensure a simplex can be built
1091 self.delaunay_triangulation(n_prc=self.n_prc)
1092 self.n_prc = self.C.shape[0]
1094 if self.disp:
1095 logging.info('Triangulation completed, evaluating all '
1096 'constraints and objective function values.')
1098 if hasattr(self, 'Tri'):
1099 self.HC.vf_to_vv(self.Tri.points, self.Tri.simplices)
1101 # Process all pools
1102 # Evaluate all constraints and functions
1103 if self.disp:
1104 logging.info('Triangulation completed, evaluating all contraints '
1105 'and objective function values.')
1107 # Evaluate all constraints and functions
1108 self.HC.V.process_pools()
1109 if self.disp:
1110 logging.info('Evaluations completed.')
1112 # feasible sampling points counted by the triangulation.py routines
1113 self.fn = self.HC.V.nfev
1114 self.n_sampled = self.nc # nevs counted in triangulation
1115 return
1117 # Hypercube minimizers
1118 def minimizers(self):
1119 """
1120 Returns the indexes of all minimizers
1121 """
1122 self.minimizer_pool = []
1123 # Note: Can implement parallelization here
1124 for x in self.HC.V.cache:
1125 in_LMC = False
1126 if len(self.LMC.xl_maps) > 0:
1127 for xlmi in self.LMC.xl_maps:
1128 if np.all(np.array(x) == np.array(xlmi)):
1129 in_LMC = True
1130 if in_LMC:
1131 continue
1133 if self.HC.V[x].minimiser():
1134 if self.disp:
1135 logging.info('=' * 60)
1136 logging.info(f'v.x = {self.HC.V[x].x_a} is minimizer')
1137 logging.info(f'v.f = {self.HC.V[x].f} is minimizer')
1138 logging.info('=' * 30)
1140 if self.HC.V[x] not in self.minimizer_pool:
1141 self.minimizer_pool.append(self.HC.V[x])
1143 if self.disp:
1144 logging.info('Neighbors:')
1145 logging.info('=' * 30)
1146 for vn in self.HC.V[x].nn:
1147 logging.info(f'x = {vn.x} || f = {vn.f}')
1149 logging.info('=' * 60)
1150 self.minimizer_pool_F = []
1151 self.X_min = []
1152 # normalized tuple in the Vertex cache
1153 self.X_min_cache = {} # Cache used in hypercube sampling
1155 for v in self.minimizer_pool:
1156 self.X_min.append(v.x_a)
1157 self.minimizer_pool_F.append(v.f)
1158 self.X_min_cache[tuple(v.x_a)] = v.x
1160 self.minimizer_pool_F = np.array(self.minimizer_pool_F)
1161 self.X_min = np.array(self.X_min)
1163 # TODO: Only do this if global mode
1164 self.sort_min_pool()
1166 return self.X_min
1168 # Local minimisation
1169 # Minimiser pool processing
1170 def minimise_pool(self, force_iter=False):
1171 """
1172 This processing method can optionally minimise only the best candidate
1173 solutions in the minimiser pool
1175 Parameters
1176 ----------
1177 force_iter : int
1178 Number of starting minimizers to process (can be sepcified
1179 globally or locally)
1181 """
1182 # Find first local minimum
1183 # NOTE: Since we always minimize this value regardless it is a waste to
1184 # build the topograph first before minimizing
1185 lres_f_min = self.minimize(self.X_min[0], ind=self.minimizer_pool[0])
1187 # Trim minimized point from current minimizer set
1188 self.trim_min_pool(0)
1190 while not self.stop_l_iter:
1191 # Global stopping criteria:
1192 self.stopping_criteria()
1194 # Note first iteration is outside loop:
1195 if force_iter:
1196 force_iter -= 1
1197 if force_iter == 0:
1198 self.stop_l_iter = True
1199 break
1201 if np.shape(self.X_min)[0] == 0:
1202 self.stop_l_iter = True
1203 break
1205 # Construct topograph from current minimizer set
1206 # (NOTE: This is a very small topograph using only the minizer pool
1207 # , it might be worth using some graph theory tools instead.
1208 self.g_topograph(lres_f_min.x, self.X_min)
1210 # Find local minimum at the miniser with the greatest Euclidean
1211 # distance from the current solution
1212 ind_xmin_l = self.Z[:, -1]
1213 lres_f_min = self.minimize(self.Ss[-1, :], self.minimizer_pool[-1])
1215 # Trim minimised point from current minimizer set
1216 self.trim_min_pool(ind_xmin_l)
1218 # Reset controls
1219 self.stop_l_iter = False
1220 return
1222 def sort_min_pool(self):
1223 # Sort to find minimum func value in min_pool
1224 self.ind_f_min = np.argsort(self.minimizer_pool_F)
1225 self.minimizer_pool = np.array(self.minimizer_pool)[self.ind_f_min]
1226 self.minimizer_pool_F = np.array(self.minimizer_pool_F)[
1227 self.ind_f_min]
1228 return
1230 def trim_min_pool(self, trim_ind):
1231 self.X_min = np.delete(self.X_min, trim_ind, axis=0)
1232 self.minimizer_pool_F = np.delete(self.minimizer_pool_F, trim_ind)
1233 self.minimizer_pool = np.delete(self.minimizer_pool, trim_ind)
1234 return
1236 def g_topograph(self, x_min, X_min):
1237 """
1238 Returns the topographical vector stemming from the specified value
1239 ``x_min`` for the current feasible set ``X_min`` with True boolean
1240 values indicating positive entries and False values indicating
1241 negative entries.
1243 """
1244 x_min = np.array([x_min])
1245 self.Y = spatial.distance.cdist(x_min, X_min, 'euclidean')
1246 # Find sorted indexes of spatial distances:
1247 self.Z = np.argsort(self.Y, axis=-1)
1249 self.Ss = X_min[self.Z][0]
1250 self.minimizer_pool = self.minimizer_pool[self.Z]
1251 self.minimizer_pool = self.minimizer_pool[0]
1252 return self.Ss
1254 # Local bound functions
1255 def construct_lcb_simplicial(self, v_min):
1256 """
1257 Construct locally (approximately) convex bounds
1259 Parameters
1260 ----------
1261 v_min : Vertex object
1262 The minimizer vertex
1264 Returns
1265 -------
1266 cbounds : list of lists
1267 List of size dimension with length-2 list of bounds for each
1268 dimension.
1270 """
1271 cbounds = [[x_b_i[0], x_b_i[1]] for x_b_i in self.bounds]
1272 # Loop over all bounds
1273 for vn in v_min.nn:
1274 for i, x_i in enumerate(vn.x_a):
1275 # Lower bound
1276 if (x_i < v_min.x_a[i]) and (x_i > cbounds[i][0]):
1277 cbounds[i][0] = x_i
1279 # Upper bound
1280 if (x_i > v_min.x_a[i]) and (x_i < cbounds[i][1]):
1281 cbounds[i][1] = x_i
1283 if self.disp:
1284 logging.info(f'cbounds found for v_min.x_a = {v_min.x_a}')
1285 logging.info(f'cbounds = {cbounds}')
1287 return cbounds
1289 def construct_lcb_delaunay(self, v_min, ind=None):
1290 """
1291 Construct locally (approximately) convex bounds
1293 Parameters
1294 ----------
1295 v_min : Vertex object
1296 The minimizer vertex
1298 Returns
1299 -------
1300 cbounds : list of lists
1301 List of size dimension with length-2 list of bounds for each
1302 dimension.
1303 """
1304 cbounds = [[x_b_i[0], x_b_i[1]] for x_b_i in self.bounds]
1306 return cbounds
1308 # Minimize a starting point locally
1309 def minimize(self, x_min, ind=None):
1310 """
1311 This function is used to calculate the local minima using the specified
1312 sampling point as a starting value.
1314 Parameters
1315 ----------
1316 x_min : vector of floats
1317 Current starting point to minimize.
1319 Returns
1320 -------
1321 lres : OptimizeResult
1322 The local optimization result represented as a `OptimizeResult`
1323 object.
1324 """
1325 # Use minima maps if vertex was already run
1326 if self.disp:
1327 logging.info(f'Vertex minimiser maps = {self.LMC.v_maps}')
1329 if self.LMC[x_min].lres is not None:
1330 logging.info(f'Found self.LMC[x_min].lres = '
1331 f'{self.LMC[x_min].lres}')
1332 return self.LMC[x_min].lres
1334 if self.callback is not None:
1335 logging.info('Callback for '
1336 'minimizer starting at {}:'.format(x_min))
1338 if self.disp:
1339 logging.info('Starting '
1340 'minimization at {}...'.format(x_min))
1342 if self.sampling_method == 'simplicial':
1343 x_min_t = tuple(x_min)
1344 # Find the normalized tuple in the Vertex cache:
1345 x_min_t_norm = self.X_min_cache[tuple(x_min_t)]
1346 x_min_t_norm = tuple(x_min_t_norm)
1347 g_bounds = self.construct_lcb_simplicial(self.HC.V[x_min_t_norm])
1348 if 'bounds' in self.min_solver_args:
1349 self.minimizer_kwargs['bounds'] = g_bounds
1350 logging.info(self.minimizer_kwargs['bounds'])
1352 else:
1353 g_bounds = self.construct_lcb_delaunay(x_min, ind=ind)
1354 if 'bounds' in self.min_solver_args:
1355 self.minimizer_kwargs['bounds'] = g_bounds
1356 logging.info(self.minimizer_kwargs['bounds'])
1358 if self.disp and 'bounds' in self.minimizer_kwargs:
1359 logging.info('bounds in kwarg:')
1360 logging.info(self.minimizer_kwargs['bounds'])
1362 # Local minimization using scipy.optimize.minimize:
1363 lres = minimize(self.func, x_min, **self.minimizer_kwargs)
1365 if self.disp:
1366 logging.info(f'lres = {lres}')
1368 # Local function evals for all minimizers
1369 self.res.nlfev += lres.nfev
1370 if 'njev' in lres:
1371 self.res.nljev += lres.njev
1372 if 'nhev' in lres:
1373 self.res.nlhev += lres.nhev
1375 try: # Needed because of the brain dead 1x1 NumPy arrays
1376 lres.fun = lres.fun[0]
1377 except (IndexError, TypeError):
1378 lres.fun
1380 # Append minima maps
1381 self.LMC[x_min]
1382 self.LMC.add_res(x_min, lres, bounds=g_bounds)
1384 return lres
1386 # Post local minimization processing
1387 def sort_result(self):
1388 """
1389 Sort results and build the global return object
1390 """
1391 # Sort results in local minima cache
1392 results = self.LMC.sort_cache_result()
1393 self.res.xl = results['xl']
1394 self.res.funl = results['funl']
1395 self.res.x = results['x']
1396 self.res.fun = results['fun']
1398 # Add local func evals to sampling func evals
1399 # Count the number of feasible vertices and add to local func evals:
1400 self.res.nfev = self.fn + self.res.nlfev
1401 return self.res
1403 # Algorithm controls
1404 def fail_routine(self, mes=("Failed to converge")):
1405 self.break_routine = True
1406 self.res.success = False
1407 self.X_min = [None]
1408 self.res.message = mes
1410 def sampled_surface(self, infty_cons_sampl=False):
1411 """
1412 Sample the function surface.
1414 There are 2 modes, if ``infty_cons_sampl`` is True then the sampled
1415 points that are generated outside the feasible domain will be
1416 assigned an ``inf`` value in accordance with SHGO rules.
1417 This guarantees convergence and usually requires less objective
1418 function evaluations at the computational costs of more Delaunay
1419 triangulation points.
1421 If ``infty_cons_sampl`` is False, then the infeasible points are
1422 discarded and only a subspace of the sampled points are used. This
1423 comes at the cost of the loss of guaranteed convergence and usually
1424 requires more objective function evaluations.
1425 """
1426 # Generate sampling points
1427 if self.disp:
1428 logging.info('Generating sampling points')
1429 self.sampling(self.nc, self.dim)
1430 if len(self.LMC.xl_maps) > 0:
1431 self.C = np.vstack((self.C, np.array(self.LMC.xl_maps)))
1432 if not infty_cons_sampl:
1433 # Find subspace of feasible points
1434 if self.g_cons is not None:
1435 self.sampling_subspace()
1437 # Sort remaining samples
1438 self.sorted_samples()
1440 # Find objective function references
1441 self.n_sampled = self.nc
1443 def sampling_custom(self, n, dim):
1444 """
1445 Generates uniform sampling points in a hypercube and scales the points
1446 to the bound limits.
1447 """
1448 # Generate sampling points.
1449 # Generate uniform sample points in [0, 1]^m \subset R^m
1450 if self.n_sampled == 0:
1451 self.C = self.sampling_function(n, dim)
1452 else:
1453 self.C = self.sampling_function(n, dim)
1454 # Distribute over bounds
1455 for i in range(len(self.bounds)):
1456 self.C[:, i] = (self.C[:, i] *
1457 (self.bounds[i][1] - self.bounds[i][0])
1458 + self.bounds[i][0])
1459 return self.C
1461 def sampling_subspace(self):
1462 """Find subspace of feasible points from g_func definition"""
1463 # Subspace of feasible points.
1464 for ind, g in enumerate(self.g_cons):
1465 # C.shape = (Z, dim) where Z is the number of sampling points to
1466 # evaluate and dim is the dimensionality of the problem.
1467 # the constraint function may not be vectorised so have to step
1468 # through each sampling point sequentially.
1469 feasible = np.array(
1470 [np.all(g(x_C, *self.g_args[ind]) >= 0.0) for x_C in self.C],
1471 dtype=bool
1472 )
1473 self.C = self.C[feasible]
1475 if self.C.size == 0:
1476 self.res.message = ('No sampling point found within the '
1477 + 'feasible set. Increasing sampling '
1478 + 'size.')
1479 # sampling correctly for both 1-D and >1-D cases
1480 if self.disp:
1481 logging.info(self.res.message)
1483 def sorted_samples(self): # Validated
1484 """Find indexes of the sorted sampling points"""
1485 self.Ind_sorted = np.argsort(self.C, axis=0)
1486 self.Xs = self.C[self.Ind_sorted]
1487 return self.Ind_sorted, self.Xs
1489 def delaunay_triangulation(self, n_prc=0):
1490 if hasattr(self, 'Tri') and self.qhull_incremental:
1491 # TODO: Uncertain if n_prc needs to add len(self.LMC.xl_maps)
1492 # in self.sampled_surface
1493 self.Tri.add_points(self.C[n_prc:, :])
1494 else:
1495 try:
1496 self.Tri = spatial.Delaunay(self.C,
1497 incremental=self.qhull_incremental,
1498 )
1499 except spatial.QhullError:
1500 if str(sys.exc_info()[1])[:6] == 'QH6239':
1501 logging.warning('QH6239 Qhull precision error detected, '
1502 'this usually occurs when no bounds are '
1503 'specified, Qhull can only run with '
1504 'handling cocircular/cospherical points'
1505 ' and in this case incremental mode is '
1506 'switched off. The performance of shgo '
1507 'will be reduced in this mode.')
1508 self.qhull_incremental = False
1509 self.Tri = spatial.Delaunay(self.C,
1510 incremental=
1511 self.qhull_incremental)
1512 else:
1513 raise
1515 return self.Tri
1518class LMap:
1519 def __init__(self, v):
1520 self.v = v
1521 self.x_l = None
1522 self.lres = None
1523 self.f_min = None
1524 self.lbounds = []
1527class LMapCache:
1528 def __init__(self):
1529 self.cache = {}
1531 # Lists for search queries
1532 self.v_maps = []
1533 self.xl_maps = []
1534 self.xl_maps_set = set()
1535 self.f_maps = []
1536 self.lbound_maps = []
1537 self.size = 0
1539 def __getitem__(self, v):
1540 try:
1541 v = np.ndarray.tolist(v)
1542 except TypeError:
1543 pass
1544 v = tuple(v)
1545 try:
1546 return self.cache[v]
1547 except KeyError:
1548 xval = LMap(v)
1549 self.cache[v] = xval
1551 return self.cache[v]
1553 def add_res(self, v, lres, bounds=None):
1554 v = np.ndarray.tolist(v)
1555 v = tuple(v)
1556 self.cache[v].x_l = lres.x
1557 self.cache[v].lres = lres
1558 self.cache[v].f_min = lres.fun
1559 self.cache[v].lbounds = bounds
1561 # Update cache size
1562 self.size += 1
1564 # Cache lists for search queries
1565 self.v_maps.append(v)
1566 self.xl_maps.append(lres.x)
1567 self.xl_maps_set.add(tuple(lres.x))
1568 self.f_maps.append(lres.fun)
1569 self.lbound_maps.append(bounds)
1571 def sort_cache_result(self):
1572 """
1573 Sort results and build the global return object
1574 """
1575 results = {}
1576 # Sort results and save
1577 self.xl_maps = np.array(self.xl_maps)
1578 self.f_maps = np.array(self.f_maps)
1580 # Sorted indexes in Func_min
1581 ind_sorted = np.argsort(self.f_maps)
1583 # Save ordered list of minima
1584 results['xl'] = self.xl_maps[ind_sorted] # Ordered x vals
1585 self.f_maps = np.array(self.f_maps)
1586 results['funl'] = self.f_maps[ind_sorted]
1587 results['funl'] = results['funl'].T
1589 # Find global of all minimizers
1590 results['x'] = self.xl_maps[ind_sorted[0]] # Save global minima
1591 results['fun'] = self.f_maps[ind_sorted[0]] # Save global fun value
1593 self.xl_maps = np.ndarray.tolist(self.xl_maps)
1594 self.f_maps = np.ndarray.tolist(self.f_maps)
1595 return results