ENH: minimize_cyipopt: support SciPy methods

cyipopt/scipy_interface.py

@@ -19,7 +19,7 @@
 else:
     SCIPY_INSTALLED = True
     del scipy
-    from scipy.optimize import approx_fprime
+    from scipy.optimize import approx_fprime, minimize
     import scipy.sparse
     try:
         from scipy.optimize import OptimizeResult
@@ -130,8 +130,10 @@ def __init__(self,
         if hess is not None:
             self.obj_hess = hess
         if jac is None:
-            jac = lambda x0, *args, **kwargs: approx_fprime(
-                x0, fun, eps, *args, **kwargs)
+            def jac(x, *args, **kwargs):
+                def wrapped_fun(x):
+                    return fun(x, *args, **kwargs)
+                return approx_fprime(x, wrapped_fun, eps)
         elif jac is True:
             fun = MemoizeJac(fun)
             jac = fun.derivative
@@ -342,6 +344,33 @@ def convert_to_bytes(options):
                 pass
 
 
+def _wrap_fun(fun, kwargs):
+    if callable(fun) and kwargs:
+        def new_fun(x, *args):
+            return fun(x, *args, **kwargs)
+    else:
+        new_fun = fun
+    return new_fun
+
+def _wrap_funs(fun, jac, hess, hessp, constraints, kwargs):
+    wrapped_fun = _wrap_fun(fun, kwargs)
+    wrapped_jac = _wrap_fun(jac, kwargs)
+    wrapped_hess = _wrap_fun(hess, kwargs)
+    wrapped_hessp = _wrap_fun(hessp, kwargs)
+    if isinstance(constraints, dict):
+        constraints = (constraints,)
+    wrapped_constraints = []
+    for constraint in constraints:
+        constraint = constraint.copy()
+        ckwargs = constraint.pop('kwargs', {})
+        constraint['fun'] = _wrap_fun(constraint.get('fun', None), ckwargs)
+        constraint['jac'] = _wrap_fun(constraint.get('jac', None), ckwargs)
+        wrapped_constraints.append(constraint)
+    return (wrapped_fun, wrapped_jac, wrapped_hess, wrapped_hessp,
+            wrapped_constraints)
+
+
+
 def minimize_ipopt(fun,
                    x0,
                    args=(),
@@ -359,6 +388,14 @@ def minimize_ipopt(fun,
     Minimization using Ipopt with an interface like
     :py:func:`scipy.optimize.minimize`.
 
+    Differences compared to :py:func:`scipy.optimize.minimize` include:
+
+    - A different default `method`: when `method` is not provided, Ipopt is
+      used to solve the problem.
+    - Support for parameter `kwargs`: additional keyword arguments to be
+      passed to the objective function, constraints, and their derivatives.
+    - Lack of support for `callback` and `hessp` with the default `method`.
+
     This function can be used to solve general nonlinear programming problems
     of the form:
 
@@ -396,8 +433,8 @@ def minimize_ipopt(fun,
         Extra keyword arguments passed to the objective function and its
         derivatives (``fun``, ``jac``, ``hess``).
     method : str, optional
-        This parameter is ignored. `minimize_ipopt` always uses Ipopt; use
-        :py:func:`scipy.optimize.minimize` directly for other methods.
+        If unspecified (default), Ipopt is used.
+        :py:func:`scipy.optimize.minimize` methods can also be used.
     jac : callable, optional
         The Jacobian of the objective function: ``jac(x, *args, **kwargs) ->
         ndarray, shape(n, )``. If ``None``, SciPy's ``approx_fprime`` is used.
@@ -406,8 +443,9 @@ def minimize_ipopt(fun,
         ``hess(x) -> ndarray, shape(n, )``.
         If ``None``, the Hessian is computed using IPOPT's numerical methods.
     hessp : callable, optional
-        This parameter is currently unused. An error will be raised if a value
-        other than ``None`` is provided.
+        If `method` is one of the SciPy methods, this is a callable that
+        produces the inner product of the Hessian and a vector. Otherwise, an
+        error will be raised if a value other than ``None`` is provided.
     bounds :  sequence, shape(n, ), optional
         Sequence of ``(min, max)`` pairs for each element in `x`. Use ``None``
         to specify no bound.
@@ -428,7 +466,8 @@ def minimize_ipopt(fun,
         and ``max_iter``. All other options are passed directly to Ipopt. See
         [1]_ for details.
     callback : callable, optional
-        This parameter is ignored.
+        This parameter is ignored unless `method` is one of the SciPy
+        methods.
 
     References
     ----------
@@ -489,6 +528,13 @@ def minimize_ipopt(fun,
     (fun, x0, args, kwargs, method, jac, hess, hessp,
      bounds, constraints, tol, callback, options) = res
 
+    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)
+        return res
+
     _x0 = np.atleast_1d(x0)
 
     lb, ub = get_bounds(bounds)
@@ -567,13 +613,14 @@ def _minimize_ipopt_iv(fun, x0, args, kwargs, method, jac, hess, hessp,
     if not np.issubdtype(x0.dtype, np.number):
         raise ValueError('`x0` must be a numeric array.')
 
-    if method is not None:  # this will be updated when gh-200 is merged
-        raise NotImplementedError('`method` is not yet supported.`')
+    # `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`)
 
-    if callback is not None:
+    if method is None and callback is not None:
         raise NotImplementedError('`callback` is not yet supported by Ipopt.`')
 
     if tol is not None:

cyipopt/tests/unit/test_scipy_optional.py

@@ -14,6 +14,11 @@
 
 import cyipopt
 
+# Hard-code rather than importing from scipy.optimize._minimize in a try/except
+try:
+    from scipy.optimize._minimize import MINIMIZE_METHODS
+except ImportError:
+    MINIMIZE_METHODS = []
 
 @pytest.mark.skipif("scipy" in sys.modules,
                     reason="Test only valid if no Scipy available.")
@@ -44,8 +49,8 @@ def f(x):
     with pytest.raises(ValueError, match=message):
         cyipopt.minimize_ipopt(f, x0, kwargs='elderberries')
 
-    message = "`method` is not yet supported."
-    with pytest.raises(NotImplementedError, match=message):
+    message = "Unknown solver..."
+    with pytest.raises(ValueError, match=message):
         cyipopt.minimize_ipopt(f, x0, method='a newt')
 
     message = "`jac` must be callable or boolean."
@@ -183,6 +188,155 @@ def test_minimize_ipopt_jac_hessians_constraints_with_arg_kwargs():
     np.testing.assert_allclose(res.get("x"), expected_res, rtol=1e-5)
 
 
+@pytest.mark.skipif("scipy" not in sys.modules,
+                    reason="Test only valid if Scipy available.")
+@pytest.mark.parametrize('method', [None] + MINIMIZE_METHODS)
+def test_minimize_ipopt_jac_with_scipy_methods(method):
+    x0 = [0] * 4
+    a0, b0, c0, d0 = 1, 2, 3, 4
+    atol, rtol = 5e-5, 5e-5
+
+    def fun(x, a=0, e=0, b=0):
+        assert a == a0
+        assert b == b0
+        fun.count += 1
+        return (x[0] - a) ** 2 + (x[1] - b) ** 2 + x[2] ** 2 + x[3] ** 2
+
+    def grad(x, a=0, e=0, b=0):
+        assert a == a0
+        assert b == b0
+        grad.count += 1
+        return [2 * (x[0] - a), 2 * (x[1] - b), 2 * x[2], 2 * x[3]]
+
+    def hess(x, a=0, e=0, b=0):
+        assert a == a0
+        assert b == b0
+        hess.count += 1
+        return 2 * np.eye(4)
+
+    def hessp(x, p, a=0, e=0, b=0):
+        assert a == a0
+        assert b == b0
+        hessp.count += 1
+        return 2 * np.eye(4) @ p
+
+    def fun_constraint(x, c=0, e=0, d=0):
+        assert c == c0
+        assert d == d0
+        fun_constraint.count += 1
+        return [(x[2] - c) ** 2, (x[3] - d) ** 2]
+
+    def grad_constraint(x, c=0, e=0, d=0):
+        assert c == c0
+        assert d == d0
+        grad_constraint.count += 1
+        return np.hstack((np.zeros((2, 2)),
+                         np.diag([2 * (x[2] - c), 2 * (x[3] - d)])))
+
+    def callback(*args, **kwargs):
+        callback.count += 1
+
+    constr = {
+        "type": "eq",
+        "fun": fun_constraint,
+        "jac": grad_constraint,
+        "args": (c0,),
+        "kwargs": {'d': d0},
+    }
+
+    kwargs = {}
+    jac_methods = {'cg', 'bfgs', 'newton-cg', 'l-bfgs-b', 'tnc', 'slsqp,',
+                   'dogleg', 'trust-ncg', 'trust-krylov', 'trust-exact',
+                   'trust-constr'}
+    hess_methods = {'newton-cg', 'dogleg', 'trust-ncg', 'trust-krylov',
+                   'trust-exact', 'trust-constr'}
+    hessp_methods = hess_methods - {'dogleg', 'trust-exact'}
+    constr_methods = {'slsqp', 'trust-constr'}
+
+    if method in jac_methods:
+        kwargs['jac'] = grad
+    if method in hess_methods:
+        kwargs['hess'] = hess
+    if method in constr_methods:
+        kwargs['constraints'] = constr
+    if method in MINIMIZE_METHODS:
+        kwargs['callback'] = callback
+
+    fun.count = 0
+    grad.count = 0
+    hess.count = 0
+    hessp.count = 0
+    fun_constraint.count = 0
+    grad_constraint.count = 0
+    callback.count = 0
+
+    res = cyipopt.minimize_ipopt(fun, x0, method=method, args=(a0,),
+                                 kwargs={'b': b0}, **kwargs)
+
+    assert res.success
+    np.testing.assert_allclose(res.x[:2], [a0, b0], rtol=rtol)
+
+    # confirm that the test covers what we think it does: all the functions
+    # that we provide are actually being executed; that is, the assertions
+    # are *passing*, not being skipped
+    assert fun.count > 0
+    if method in MINIMIZE_METHODS:
+        assert callback.count > 0
+    if method in jac_methods:
+        assert grad.count > 0
+    if method in hess_methods:
+        assert hess.count > 0
+    if method in constr_methods:
+        assert fun_constraint.count > 0
+        assert grad_constraint.count > 0
+        np.testing.assert_allclose(res.x[2:], [c0, d0], rtol=rtol)
+    else:
+        np.testing.assert_allclose(res.x[2:], 0, atol=atol)
+
+    # For methods that support `hessp`, check that it works, too.
+    if method in hessp_methods:
+        del kwargs['hess']
+        kwargs['hessp'] = hessp
+
+        res = cyipopt.minimize_ipopt(fun, x0, method=method, args=(a0,),
+                                     kwargs={'b': b0}, **kwargs)
+        assert res.success
+        assert hessp.count > 0
+        np.testing.assert_allclose(res.x[:2], [a0, b0], rtol=rtol)
+        if method in constr_methods:
+            np.testing.assert_allclose(res.x[2:], [c0, d0], rtol=rtol)
+        else:
+            np.testing.assert_allclose(res.x[2:], 0, atol=atol)
+
+
+@pytest.mark.skipif("scipy" not in sys.modules,
+                    reason="Test only valid of Scipy available")
+def test_minimize_ipopt_bounds_tol_options():
+    # spot check additional cases not tested above
+    def fun(x):
+        return x**2
+
+    x0 = 2.
+
+    # make sure `bounds` is passed to SciPy methods
+    bounds = (1, 3)
+    res = cyipopt.minimize_ipopt(fun, x0, method='slsqp', bounds=[(1, 3)])
+    np.testing.assert_allclose(res.x, bounds[0])
+
+    # make sure `tol` is passed to SciPy methods
+    tol1, tol2 = 1e-3, 1e-9
+    res1 = cyipopt.minimize_ipopt(fun, x0, method='cobyla', tol=tol1)
+    res2 = cyipopt.minimize_ipopt(fun, x0, method='cobyla', tol=tol2)
+    assert abs(res2.x[0]) <= tol2 < abs(res1.x[0]) <= tol1
+
+    # make sure `options` is passed to SciPy methods
+    res = cyipopt.minimize_ipopt(fun, x0, method='slsqp')
+    assert res.nit > 1
+    options = dict(maxiter=1)
+    res = cyipopt.minimize_ipopt(fun, x0, method='slsqp', options=options)
+    assert res.nit == 1
+
+
 @pytest.mark.skipif("scipy" not in sys.modules,
                     reason="Test only valid of Scipy available")
 def test_minimize_ipopt_sparse_jac_if_scipy():