https://github.com/lmfit/lmfit-py
Tip revision: 9f9af6f36c0928767ea8b004ea8cb5a16aba6b04 authored by Matt Newville on 14 October 2021, 19:34:30 UTC
Merge pull request #759 from reneeotten/last_changes_before_release
Merge pull request #759 from reneeotten/last_changes_before_release
Tip revision: 9f9af6f
_ampgo.py
"""Adaptive Memory Programming for Global Optimization (AMPGO).
added to lmfit by Renee Otten (2018)
based on the Python implementation of Andrea Gavana
(see: http://infinity77.net/global_optimization/)
Implementation details can be found in this paper:
http://leeds-faculty.colorado.edu/glover/fred%20pubs/416%20-%20AMP%20(TS)%20for%20Constrained%20Global%20Opt%20w%20Lasdon%20et%20al%20.pdf
"""
import numpy as np
from scipy.optimize import minimize
SCIPY_LOCAL_SOLVERS = ['Nelder-Mead', 'Powell', 'L-BFGS-B', 'TNC', 'SLSQP']
def ampgo(objfun, x0, args=(), local='L-BFGS-B', local_opts=None, bounds=None,
maxfunevals=None, totaliter=20, maxiter=5, glbtol=1e-5, eps1=0.02,
eps2=0.1, tabulistsize=5, tabustrategy='farthest', disp=False):
"""Find the global minimum of a multivariate function using AMPGO.
AMPGO stands for "Adaptive Memory Programming for Global Optimization.
Parameters
----------
objfun : callable
Objective function to be minimized. The function must have the
signature::
objfun(params, *args, **kws)
x0 : numpy.ndarray
Initial guesses for parameter values.
args : tuple, optional
Additional arguments passed to `objfun`.
local : str, optional
Name of the local minimization method. Valid options are:
- `'L-BFGS-B'` (default)
- `'Nelder-Mead'`
- `'Powell'`
- `'TNC'`
- `'SLSQP'`
local_opts : dict, optional
Options to pass to the local minimizer.
bounds : sequence, optional
List of tuples specifying the lower and upper bound for each
independent variable ``[(xl0, xu0), (xl1, xu1), ...]``.
maxfunevals : int, optional
Maximum number of function evaluations. If None, the optimization
will stop after `totaliter` number of iterations.
totaliter : int, optional
Maximum number of global iterations.
maxiter : int, optional
Maximum number of 'Tabu Tunneling' iterations during each global
iteration.
glbtol : float, optional
Tolerance whether or not to accept a solution after a tunneling
phase.
eps1 : float, optional
Constant used to define an aspiration value for the objective
function during the Tunneling phase.
eps2 : float, optional
Perturbation factor used to move away from the latest local
minimum at the start of a Tunneling phase.
tabulistsize : int, optional
Size of the (circular) tabu search list.
tabustrategy : {'farthest', 'oldest'}, optional
Strategy to use when the size of the tabu list exceeds
`tabulistsize`. It can be 'oldest' to drop the oldest point from
the tabu list or 'farthest' (default) to drop the element farthest
from the last local minimum found.
disp : bool, optional
Set to True to print convergence messages.
Returns
-------
tuple
A tuple of 5 elements, in the following order:
1. **best_x** (array_like): the estimated position of the global
minimum.
2. **best_f** (float): the value of `objfun` at the minimum.
3. **evaluations** (int): the number of function evaluations.
4. **msg** (str): a message describes the cause of the
termination.
5. **tunnel_info** (tuple): a tuple containing the total number
of Tunneling phases performed and the successful ones.
Notes
-----
The detailed implementation of AMPGO is described in the paper
"Adaptive Memory Programming for Constrained Global Optimization"
located here:
http://leeds-faculty.colorado.edu/glover/fred%20pubs/416%20-%20AMP%20(TS)%20for%20Constrained%20Global%20Opt%20w%20Lasdon%20et%20al%20.pdf
"""
if local not in SCIPY_LOCAL_SOLVERS:
raise Exception(f'Invalid local solver selected: {local}')
x0 = np.atleast_1d(x0)
n = len(x0)
if bounds is None:
bounds = [(None, None)] * n
if len(bounds) != n:
raise ValueError('length of x0 != length of bounds')
bounds = [b if b is not None else (None, None) for b in bounds]
_bounds = [(-np.inf if lb is None else lb, np.inf if ub is None else ub)
for lb, ub in bounds]
low, up = zip(*_bounds)
if maxfunevals is None:
maxfunevals = np.inf
if tabulistsize < 1:
raise Exception(f'Invalid tabulistsize specified: {tabulistsize}. '
'It should be an integer greater than zero.')
if tabustrategy not in ['oldest', 'farthest']:
raise Exception(f'Invalid tabustrategy specified: {tabustrategy}. '
'It must be one of "oldest" or "farthest".')
tabulist = []
best_f = np.inf
best_x = x0
global_iter = 0
all_tunnel = success_tunnel = 0
evaluations = 0
local_tol = min(1e-8, glbtol)
while 1:
# minimization to find local minimum, either from initial values or
# after a successful tunneling loop
if disp:
print('\n{0}\nStarting MINIMIZATION Phase {1:d}\n{0}'
.format('='*72, global_iter+1))
options = {'maxiter': max(1, maxfunevals), 'disp': disp}
if local_opts is not None:
options.update(local_opts)
res = minimize(objfun, x0, args=args, method=local, bounds=bounds,
tol=local_tol, options=options)
xf, yf, num_fun = res['x'], res['fun'], res['nfev']
if isinstance(yf, np.ndarray):
yf = yf[0]
maxfunevals -= num_fun
evaluations += num_fun
if yf < best_f:
best_f = yf
best_x = xf
if disp:
print(f'\n\n ==> Reached local minimum: {yf:.5g}\n')
if maxfunevals <= 0:
if disp:
print('='*72)
return (best_x, best_f, evaluations,
'Maximum number of function evaluations exceeded',
(all_tunnel, success_tunnel))
# if needed, drop a value from the tabu tunneling list and add the
# current solution
tabulist = drop_tabu_points(xf, tabulist, tabulistsize, tabustrategy)
tabulist.append(xf)
i = improve = 0
while i < maxiter and improve == 0:
if disp:
print('{0}\nStarting TUNNELING Phase ({1:d}-{2:d})\n{0}'
.format('='*72, global_iter+1, i+1))
all_tunnel += 1
# generate a new starting point away from the current solution
r = np.random.uniform(-1.0, 1.0, size=(n, ))
beta = eps2*np.linalg.norm(xf) / np.linalg.norm(r)
if np.abs(beta) < 1e-8:
beta = eps2
x0 = xf + beta*r
# make sure that the new starting point is within bounds
x0 = np.where(x0 < low, low, x0)
x0 = np.where(x0 > up, up, x0)
# aspired value of the objective function for the tunneling loop
aspiration = best_f - eps1*(1.0 + np.abs(best_f))
tunnel_args = tuple([objfun, aspiration, tabulist] + list(args))
options = {'maxiter': max(1, maxfunevals), 'disp': disp}
if local_opts is not None:
options.update(local_opts)
res = minimize(tunnel, x0, args=tunnel_args, method=local,
bounds=bounds, tol=local_tol, options=options)
xf, yf, num_fun = res['x'], res['fun'], res['nfev']
if isinstance(yf, np.ndarray):
yf = yf[0]
maxfunevals -= num_fun
evaluations += num_fun
yf = inverse_tunnel(xf, yf, aspiration, tabulist)
if yf <= best_f + glbtol:
oldf = best_f
best_f = yf
best_x = xf
improve = 1
success_tunnel += 1
if disp:
print('\n\n ==> Successful tunnelling phase. Reached new '
f'local minimum: {yf:.5g} < {oldf:.5g}\n')
i += 1
if maxfunevals <= 0:
return (best_x, best_f, evaluations,
'Maximum number of function evaluations exceeded',
(all_tunnel, success_tunnel))
tabulist = drop_tabu_points(xf, tabulist, tabulistsize, tabustrategy)
tabulist.append(xf)
if disp:
print('='*72)
global_iter += 1
x0 = xf.copy()
if global_iter >= totaliter:
return (best_x, best_f, evaluations,
'Maximum number of global iterations exceeded',
(all_tunnel, success_tunnel))
def drop_tabu_points(xf, tabulist, tabulistsize, tabustrategy):
"""Drop a point from the tabu search list."""
if len(tabulist) < tabulistsize:
return tabulist
if tabustrategy == 'oldest':
tabulist.pop(0)
else:
distance = np.sqrt(np.sum((tabulist - xf)**2, axis=1))
index = np.argmax(distance)
tabulist.pop(index)
return tabulist
def tunnel(x0, *args):
"""Tunneling objective function.
This function has a global minimum of zero at any feasible point where
``f(x) = aspiration``, and minimizing this expression tends to move
away from all points in `tabulist`.
"""
objfun, aspiration, tabulist, *fun_args = args
numerator = (objfun(x0, *fun_args) - aspiration)**2
denominator = 1.0
for tabu in tabulist:
denominator = denominator * np.sqrt(np.sum((x0 - tabu)**2))
ytf = numerator/denominator
return ytf
def inverse_tunnel(xtf, ytf, aspiration, tabulist):
"""Calculate the function value after a tunneling phase step."""
denominator = 1.0
for tabu in tabulist:
denominator = denominator * np.sqrt(np.sum((xtf - tabu)**2))
numerator = ytf*denominator
yf = aspiration + np.sqrt(numerator)
return yf