262 lines
9.4 KiB
Python
262 lines
9.4 KiB
Python
from .utils import wrap_untrustworthy, check, final, ExhaustedTrialsError
|
|
import numpy as np
|
|
import scipy.optimize as scopt
|
|
|
|
|
|
def make_scipy(method, *, jacobian=None, hessian=None):
|
|
# does not include:
|
|
# basinhopping, shgo
|
|
assert method in (
|
|
"Nelder-Mead",
|
|
"Powell",
|
|
"CG",
|
|
"BFGS",
|
|
# "Newton-CG",
|
|
"L-BFGS-B",
|
|
"TNC",
|
|
"COBYLA",
|
|
"SLSQP",
|
|
"trust-constr",
|
|
# "dogleg",
|
|
# "trust-ncg",
|
|
# "trust-exact",
|
|
# "trust-krylov",
|
|
), method
|
|
|
|
def f(objective, n_trials, n_dim, with_count):
|
|
prng = np.random.default_rng()
|
|
_objective = wrap_untrustworthy(
|
|
objective, n_trials, raising=True, bounding="sine"
|
|
)
|
|
x0 = np.full(n_dim, 0.5)
|
|
bounds = scopt.Bounds([0.0] * n_dim, [1.0] * n_dim)
|
|
# jac = "cs" if method == "dogleg" else None # doesn't work
|
|
jac = "2-point" if jacobian is None else jacobian
|
|
hess = "2-point" if hessian is None else hessian
|
|
# Alternatively, objects implementing the HessianUpdateStrategy interface
|
|
# can be used to approximate the Hessian. Available quasi-Newton methods
|
|
# implementing this interface are: BFGS; SR1.
|
|
tol = None # 0.0
|
|
# options = dict(maxfun=n_trials) if method == "TNC" else dict(maxiter=n_trials)
|
|
|
|
if method in ("BFGS", "CG", "COBYLA", "SLSQP", "trust-constr"):
|
|
options = dict(maxiter=n_trials)
|
|
elif method in ("Nelder-Mead", "Powell"):
|
|
options = dict(maxfev=n_trials, maxiter=n_trials)
|
|
elif method in ("L-BFGS-B", "TNC"):
|
|
options = dict(maxfun=n_trials, maxiter=n_trials)
|
|
else:
|
|
options = dict(maxfun=n_trials, maxfev=n_trials, maxiter=n_trials)
|
|
|
|
# silence some warnings:
|
|
if method in ("Nelder-Mead", "Powell"):
|
|
jac = None
|
|
hess = None
|
|
if method in ("COBYLA",):
|
|
jac = None
|
|
hess = None
|
|
bounds = None
|
|
if method in ("BFGS", "CG"):
|
|
hess = None
|
|
bounds = None
|
|
if method in ("L-BFGS-B", "SLSQP", "TNC", "trust-constr"):
|
|
hess = None
|
|
|
|
checks = []
|
|
|
|
def check_evals():
|
|
evals = _objective(check)
|
|
checks.append(evals)
|
|
return evals < n_trials
|
|
|
|
first_try = True
|
|
while check_evals():
|
|
if not first_try:
|
|
x0 = prng.uniform(size=n_dim)
|
|
|
|
try:
|
|
res = scopt.minimize(
|
|
_objective,
|
|
x0,
|
|
method=method,
|
|
jac=jac,
|
|
hess=hess,
|
|
bounds=bounds,
|
|
tol=tol,
|
|
options=options,
|
|
)
|
|
except ExhaustedTrialsError:
|
|
break
|
|
else:
|
|
# well, this is pointless.
|
|
fopt, xopt, feval_count = res.fun, res.x, res.nfev
|
|
# print("success:", res.success)
|
|
# if not res.success and res.nfev < n_trials // 2:
|
|
shut_up = (
|
|
method == "SLSQP"
|
|
and res.message == "Inequality constraints incompatible"
|
|
) or (
|
|
res.message
|
|
== "The maximum number of function evaluations is exceeded."
|
|
)
|
|
if not shut_up and not res.success and res.nfev < n_trials // 3:
|
|
print("", method, res, "", sep="\n")
|
|
|
|
first_try = False
|
|
|
|
# TODO: run without this, try to minimize number of attempts (i.e. list length)
|
|
# if len(checks) >= 5: print(method, [b - a for a, b in zip(checks, checks[1:])])
|
|
|
|
fopt, xopt, feval_count = _objective(final)
|
|
return (fopt, xopt, feval_count) if with_count else (fopt, xopt)
|
|
|
|
name = f"scipy_{method.replace('-', '').lower()}"
|
|
if jacobian == "2-point":
|
|
name += "_2j"
|
|
elif jacobian == "3-point":
|
|
name += "_3j"
|
|
elif jacobian is not None:
|
|
assert False, jacobian
|
|
if hessian == "2-point":
|
|
name += "_2h"
|
|
elif hessian == "3-point":
|
|
name += "_3h"
|
|
elif hessian is not None:
|
|
assert False, hessian
|
|
f.__name__ = name + "_cube"
|
|
return f
|
|
|
|
|
|
def scipy_basinhopping_cube(objective, n_trials, n_dim, with_count):
|
|
progress = 1e-2 # TODO: make configurable?
|
|
# NOTE: could also callbacks to extract solutions instead of wrapping objective functions?
|
|
|
|
def accept_bounded(x_new=None, x_old=None, f_new=None, f_old=None):
|
|
return np.all(x_new >= 0.0) and np.all(x_new <= 1.0)
|
|
|
|
def dummy_minimizer(fun, x0, args, **options):
|
|
return scopt.OptimizeResult(x=x0, fun=fun(x0), success=True, nfev=1)
|
|
|
|
x0 = np.full(n_dim, 0.5)
|
|
res = scopt.basinhopping(
|
|
objective,
|
|
x0,
|
|
minimizer_kwargs=dict(method=dummy_minimizer),
|
|
accept_test=accept_bounded,
|
|
disp=False,
|
|
niter=n_trials,
|
|
# TODO: try without any progress vars at all.
|
|
T=progress,
|
|
stepsize=progress / 2,
|
|
)
|
|
fopt, xopt, feval_count = res.fun, res.x, res.nfev
|
|
# print("success:", res.success)
|
|
return (fopt, xopt, feval_count) if with_count else (fopt, xopt)
|
|
|
|
|
|
def scipy_direct_cube(objective, n_trials, n_dim, with_count):
|
|
bounds = scopt.Bounds([0.0] * n_dim, [1.0] * n_dim)
|
|
# TODO: try different values of eps. default 0.0001
|
|
res = scopt.direct(
|
|
objective,
|
|
bounds=bounds,
|
|
maxfun=n_trials,
|
|
maxiter=1_000_000,
|
|
vol_tol=0,
|
|
)
|
|
fopt, xopt, feval_count = res.fun, res.x, res.nfev
|
|
# print("success:", res.success)
|
|
return (fopt, xopt, feval_count) if with_count else (fopt, xopt)
|
|
|
|
|
|
def scipy_direct_l_cube(objective, n_trials, n_dim, with_count):
|
|
bounds = scopt.Bounds([0.0] * n_dim, [1.0] * n_dim)
|
|
# TODO: try different values of eps. default 0.0001
|
|
res = scopt.direct(
|
|
objective,
|
|
bounds=bounds,
|
|
maxfun=n_trials,
|
|
maxiter=1_000_000,
|
|
vol_tol=0,
|
|
locally_biased=True,
|
|
len_tol=0.0, # only for locally_biased=True
|
|
)
|
|
fopt, xopt, feval_count = res.fun, res.x, res.nfev
|
|
# print("success:", res.success)
|
|
return (fopt, xopt, feval_count) if with_count else (fopt, xopt)
|
|
|
|
|
|
def make_shgo(method="cobyla", init="large", it=1, mei=False, li=False, ftol=12):
|
|
from scipy.optimize import shgo, Bounds
|
|
|
|
if method == "cobyqa":
|
|
from cobyqa import minimize as cobyqa
|
|
from cobyqa.optimize import EXIT_RHOEND_SUCCESS, EXIT_MAXFEV_WARNING
|
|
|
|
raise NotImplementedError("TODO")
|
|
else:
|
|
# from scipy.optimize import minimize
|
|
# https://docs.scipy.org/doc/scipy/reference/optimize.minimize-cobyla.html
|
|
pass
|
|
|
|
def f(objective, n_trials, n_dim, with_count):
|
|
_objective = wrap_untrustworthy(
|
|
objective, n_trials, bounding="clip", raising=True
|
|
)
|
|
bounds = Bounds([0.0] * n_dim, [1.0] * n_dim)
|
|
npt = dict(
|
|
large=(n_dim + 1) * (n_dim + 2) // 2, medium=2 * n_dim + 1, small=n_dim + 1
|
|
)[init]
|
|
_method = method if method != "cobyqa" else method # TODO: handle cobyqa case.
|
|
try:
|
|
shgo(
|
|
_objective,
|
|
bounds=bounds,
|
|
n=npt,
|
|
iters=it,
|
|
minimizer_kwargs=dict(method=_method, ftol=10**-ftol),
|
|
options=dict(maxfev=n_trials, minimize_every_iter=mei, local_iters=li),
|
|
sampling_method="simplicial",
|
|
)
|
|
except ExhaustedTrialsError:
|
|
pass
|
|
fopt, xopt, feval_count = _objective(final)
|
|
return (fopt, xopt, feval_count) if with_count else (fopt, xopt)
|
|
|
|
name = f"shgo_{method}"
|
|
if it != 1:
|
|
name += f"_it{it}"
|
|
if li:
|
|
name += f"_li{li}"
|
|
if mei:
|
|
name += "_mei"
|
|
if init != "large":
|
|
name += f"_{init}"
|
|
if ftol != 12:
|
|
name += f"_ft{ftol:02}"
|
|
f.__name__ = name + "_cube"
|
|
return f
|
|
|
|
|
|
scipy_bfgs_2j_cube = make_scipy("BFGS", jacobian="2-point")
|
|
scipy_bfgs_3j_cube = make_scipy("BFGS", jacobian="3-point")
|
|
scipy_cg_2j_cube = make_scipy("CG", jacobian="2-point")
|
|
scipy_cg_3j_cube = make_scipy("CG", jacobian="3-point")
|
|
scipy_cobyla_cube = make_scipy("COBYLA")
|
|
# scipy_dogleg_cube = make_scipy("dogleg") # ValueError: Jacobian is required for dogleg minimization
|
|
scipy_lbfgsb_2j_cube = make_scipy("L-BFGS-B", jacobian="2-point")
|
|
scipy_lbfgsb_3j_cube = make_scipy("L-BFGS-B", jacobian="3-point")
|
|
scipy_neldermead_cube = make_scipy("Nelder-Mead")
|
|
# scipy_newtoncg_cube = make_scipy("Newton-CG") # ValueError: Jacobian is required for Newton-CG method
|
|
scipy_powell_cube = make_scipy("Powell")
|
|
scipy_slsqp_2j_cube = make_scipy("SLSQP", jacobian="2-point")
|
|
scipy_slsqp_3j_cube = make_scipy("SLSQP", jacobian="3-point")
|
|
scipy_tnc_2j_cube = make_scipy("TNC", jacobian="2-point")
|
|
scipy_tnc_3j_cube = make_scipy("TNC", jacobian="3-point")
|
|
scipy_trustconstr_2j_cube = make_scipy("trust-constr", jacobian="2-point")
|
|
scipy_trustconstr_3j_cube = make_scipy("trust-constr", jacobian="3-point")
|
|
# scipy_trustexact_2j_cube = make_scipy("trust-exact", jacobian="2-point")
|
|
# scipy_trustexact_3j_cube = make_scipy("trust-exact", jacobian="3-point")
|
|
# scipy_trustkrylov_cube = make_scipy("trust-krylov") # ValueError: ('Jacobian is required for trust region ', 'exact minimization.')
|
|
# scipy_trustncg_cube = make_scipy("trust-ncg") # ValueError: Jacobian is required for Newton-CG trust-region minimization
|