Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

ENH: minimize_ipopt: add input validation for bounds and support for new-style constraints #207

Merged
merged 14 commits into from
Sep 19, 2023
Merged
Show file tree
Hide file tree
Changes from 11 commits
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
71 changes: 53 additions & 18 deletions cyipopt/scipy_interface.py
Original file line number Diff line number Diff line change
Expand Up @@ -19,7 +19,7 @@
else:
SCIPY_INSTALLED = True
del scipy
from scipy.optimize import approx_fprime, minimize
from scipy import optimize
Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I needed some other things from optimize.

import scipy.sparse
try:
from scipy.optimize import OptimizeResult
Expand Down Expand Up @@ -129,11 +129,11 @@ def __init__(self,

if hess is not None:
self.obj_hess = hess
if jac is None:
if not jac:
Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

When jac is False, the code currently treats it as though it were a callable.

def jac(x, *args, **kwargs):
def wrapped_fun(x):
return fun(x, *args, **kwargs)
return approx_fprime(x, wrapped_fun, eps)
return optimize.approx_fprime(x, wrapped_fun, eps)
elif jac is True:
fun = MemoizeJac(fun)
jac = fun.derivative
Expand All @@ -159,11 +159,13 @@ def wrapped_fun(x):
con_hessian = con.get('hess', None)
con_kwargs = con.get('kwargs', {})
if con_jac is None:
con_jac = lambda x0, *args, **kwargs: optimize.approx_fprime(
x0, con_fun, eps, *args, **kwargs)
Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Suggested change
con_jac = lambda x0, *args, **kwargs: optimize.approx_fprime(
x0, con_fun, eps, *args, **kwargs)

This looks like it snuck in with a merge or something.

# beware of late binding!
def con_jac(x, *args, con_fun=con_fun, **kwargs):
def wrapped(x):
return con_fun(x, *args, **kwargs)
return approx_fprime(x, wrapped, eps)
return optimize.approx_fprime(x, wrapped, eps)
elif con_jac is True:
con_fun = MemoizeJac(con_fun)
con_jac = con_fun.derivative
Expand Down Expand Up @@ -531,17 +533,18 @@ def minimize_ipopt(fun,
if method is not None:
funs = _wrap_funs(fun, jac, hess, hessp, constraints, kwargs)
fun, jac, hess, hessp, constraints = funs
res = minimize(fun, x0, args, method, jac, hess, hessp,
bounds, constraints, tol, callback, options)
bounds = optimize.Bounds(*bounds)
res = optimize.minimize(fun, x0, args, method, jac, hess, hessp,
bounds, constraints, tol, callback, options)
return res

_x0 = np.atleast_1d(x0)
Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Conversion to 1d is now done in _minimize_ipopt_iv.

Previously, the variable was renamed to _x0 so that at the end, if x0 had been a scalar, the x attribute of the result object would be a scalar.

In SciPy, I am very careful about shapes, and I have corrected many shape bugs related to this sort of thing in SciPy code. However, I'm not sure if this is needed here.

The documentation (which now mirrors that of scipy.optimize.minimize) states that x0 should be an array-like of shape (n,), where n is the number of independent variables. This function does not handle arrays of arbitrary dimensionality, and there is a separate minimize_scalar function for functions of a single variable. Also, the x attribute of the object returned by scipy.optimize.minimize returns 1d array even when x0 is a scalar:

from scipy.optimize import minimize

def f(x):
    return x**2

minimize(f, 1).x  # [-7.450e-09]

All these things suggest to me that the output should always be a 1d array (even if the function is flexible about accepting scalar input rather than raising an error)and working with functions that return scalars).

LMK if you'd prefer for this to be left alone.

Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

If our function here supported scalar and array objective functions simultaneously, I think we'd need to deprecate the scalar functionality with at least a release.

Copy link
Contributor Author

@mdhaber mdhaber Sep 15, 2023

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

This would not change what types of objective functions are supported. It makes the function more consistent with the documentation, which states:

image

In SciPy, the x attribute of the result object is always 1D. However, currently in minimize_ipopt, the x attribute of the result object may be a scalar. The proposed change would resolve this inconsistency and make minimize_ipopt behave as advertised - like scipy.optimize.minimize.

It may surprise users who relied on the result attribute x being a scalar in some cases, but this behavior wasn't consistent with the documentation, so it could be considered a bug fix.

Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Ok. I consider this backwards incompatible, regardless of what the docs say it should do, but we can make the change with no deprecation.


lb, ub = get_bounds(bounds)
lb, ub = bounds
cl, cu = get_constraint_bounds(constraints, _x0)
con_dims = get_constraint_dimensions(constraints, _x0)
sparse_jacs, jac_nnz_row, jac_nnz_col = _get_sparse_jacobian_structure(
constraints, _x0)
constraints, x0)

problem = IpoptProblemWrapper(fun,
args=args,
Expand All @@ -559,7 +562,7 @@ def minimize_ipopt(fun,
if options is None:
options = {}

nlp = cyipopt.Problem(n=len(_x0),
nlp = cyipopt.Problem(n=len(x0),
m=len(cl),
problem_obj=problem,
lb=lb,
Expand All @@ -573,7 +576,9 @@ def minimize_ipopt(fun,
# Rename some default scipy options
replace_option(options, b'disp', b'print_level')
replace_option(options, b'maxiter', b'max_iter')
if b'print_level' not in options:
if getattr(options, 'print_level', False) is True:
Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

When disp is True, we get an error because print_level needs to be a bool.

options[b'print_level'] = 1
else:
options[b'print_level'] = 0
if b'tol' not in options:
options[b'tol'] = tol or 1e-8
Expand All @@ -589,10 +594,7 @@ def minimize_ipopt(fun,
msg = 'Invalid option for IPOPT: {0}: {1} (Original message: "{2}")'
raise TypeError(msg.format(option, value, e))

x, info = nlp.solve(_x0)

if np.asarray(x0).shape == ():
x = x[0]
x, info = nlp.solve(x0)

return OptimizeResult(x=x,
success=info['status'] == 0,
Expand All @@ -604,21 +606,54 @@ def minimize_ipopt(fun,
njev=problem.njev,
nit=problem.nit)


def _minimize_ipopt_iv(fun, x0, args, kwargs, method, jac, hess, hessp,
bounds, constraints, tol, callback, options):
# basic input validation for minimize_ipopt that is not included in
# IpoptProblemWrapper

x0 = np.asarray(x0)[()]
x0 = np.atleast_1d(x0)
if not np.issubdtype(x0.dtype, np.number):
raise ValueError('`x0` must be a numeric array.')

# `method` does not need input validation. If `method is not None`, it is
# passed to `scipy.optimize.minimize`, which raises a readable error if
# the value isn't recognized.

# TODO: add input validation for `bounds` when adding
# support for instances of new-style constraints (e.g. `Bounds`)
# Handle bounds that are either sequences (of sequences) or instances of
mdhaber marked this conversation as resolved.
Show resolved Hide resolved
# `optimize.Bounds`.
if bounds is None:
bounds = [-np.inf, np.inf]

if isinstance(bounds, optimize.Bounds):
lb, ub = bounds.lb, bounds.ub
else:
bounds = np.atleast_2d(bounds)
if bounds.shape[1] != 2:
raise ValueError("`bounds` must specify both lower and upper "
"limits for each decision variable.")
lb, ub = bounds.T

try:
lb, ub, x0 = np.broadcast_arrays(lb, ub, x0)
mdhaber marked this conversation as resolved.
Show resolved Hide resolved
except ValueError:
raise ValueError("The number of lower bounds, upper bounds, and "
"decision variables must be equal or broadcastable.")

try:
lb = lb.astype(np.float64)
ub = ub.astype(np.float64)
mdhaber marked this conversation as resolved.
Show resolved Hide resolved
except ValueError:
raise ValueError("The bounds must be numeric.")

# Nones turn into NaNs above. Previously, NaNs caused Ipopt to hang, so
# I'm not concerned about turning them into infs.
lb[np.isnan(lb)] = -np.inf
ub[np.isnan(ub)] = np.inf

bounds = lb, ub

constraints = optimize._minimize.standardize_constraints(constraints, x0,
'old')

if method is None and callback is not None:
raise NotImplementedError('`callback` is not yet supported by Ipopt.`')
Expand Down
Loading