Add CyIpoptEvaluationError and handle in callbacks

cyipopt/__init__.py

@@ -13,3 +13,4 @@
 from .ipopt_wrapper import *
 from .scipy_interface import *
 from .version import __version__
+from .exceptions import *

cyipopt/cython/ipopt_wrapper.pyx

@@ -17,6 +17,7 @@ import inspect
 import numpy as np
 cimport numpy as np
 
+from cyipopt.exceptions import CyIpoptEvaluationError
 from cyipopt.utils import deprecated_warning, generate_deprecation_warning_msg
 from ipopt cimport *
 
@@ -868,6 +869,8 @@ cdef Bool objective_cb(Index n,
         _x[i] = x[i]
     try:
         obj_value[0] = self.__objective(_x)
+    except CyIpoptEvaluationError:
+        return False
     except:
         self.__exception = sys.exc_info()
     return True
@@ -892,6 +895,8 @@ cdef Bool gradient_cb(Index n,
 
     try:
         ret_val = self.__gradient(_x)
+    except CyIpoptEvaluationError:
+        return False
     except:
         self.__exception = sys.exc_info()
         return True
@@ -928,6 +933,8 @@ cdef Bool constraints_cb(Index n,
 
     try:
         ret_val = self.__constraints(_x)
+    except CyIpoptEvaluationError:
+        return False
     except:
         self.__exception = sys.exc_info()
         return True
@@ -1002,6 +1009,8 @@ cdef Bool jacobian_cb(Index n,
 
         try:
             ret_val = self.__jacobian(_x)
+        except CyIpoptEvaluationError:
+            return False
         except:
             self.__exception = sys.exc_info()
             return True
@@ -1087,6 +1096,8 @@ cdef Bool hessian_cb(Index n,
 
         try:
             ret_val = self.__hessian(_x, _lambda, obj_factor)
+        except CyIpoptEvaluationError:
+            return False
         except:
             self.__exception = sys.exc_info()
             return True

cyipopt/exceptions.py

@@ -0,0 +1,45 @@
+# -*- coding: utf-8 -*-
+"""
+cyipopt: Python wrapper for the Ipopt optimization package, written in Cython.
+
+Copyright (C) 2012-2015 Amit Aides
+Copyright (C) 2015-2017 Matthias Kümmerer
+Copyright (C) 2017-2023 cyipopt developers
+
+License: EPL 2.0
+"""
+
+class CyIpoptEvaluationError(ArithmeticError):
+    """An exception that should be raised in evaluation callbacks to signal
+    to CyIpopt that a numerical error occured during function evaluation.
+
+    Whereas most exceptions that occur in callbacks are re-raised, exceptions
+    of this type are ignored other than to communicate to Ipopt that an error
+    occurred.
+
+    Ipopt handles evaluation errors differently depending on where they are
+    raised (which evaluation callback returns ``false`` to Ipopt).
+    When evaluation errors are raised in the following callbacks, Ipopt
+    attempts to recover by cutting the step size. This is usually the desired
+    behavior when an undefined value is encountered.
+
+    - ``objective``
+    - ``constraints``
+
+    When raised in the following callbacks, Ipopt fails with an "Invalid number"
+    return status.
+
+    - ``gradient``
+    - ``jacobian``
+    - ``hessian``
+
+    Raising an evaluation error in the following callbacks results is not
+    supported.
+
+    - ``jacobianstructure``
+    - ``hessianstructure``
+    - ``intermediate``
+
+    """
+
+    pass

cyipopt/tests/integration/test_hs071.py

@@ -15,3 +15,243 @@ def test_hs071_solve(hs071_initial_guess_fixture, hs071_problem_instance_fixture
 
 	expected_x = np.array([1.0, 4.74299964, 3.82114998, 1.37940829])
 	np.testing.assert_allclose(x, expected_x)
+
+
+def _make_problem(definition, lb, ub, cl, cu):
+    n = len(lb)
+    m = len(cl)
+    return cyipopt.Problem(
+        n=n, m=m, problem_obj=definition, lb=lb, ub=ub, cl=cl, cu=cu
+    )
+
+
+def _solve_and_assert_correct(problem, x0):
+    x, info = problem.solve(x0)
+    expected_x = np.array([1.0, 4.74299964, 3.82114998, 1.37940829])
+    assert info["status"] == 0
+    np.testing.assert_allclose(x, expected_x)
+
+
+def _assert_solve_fails(problem, x0):
+    x, info = problem.solve(x0)
+    # The current (Ipopt 3.14.11) return status is "Invalid number".
+    assert info["status"] < 0
+
+
+def test_hs071_objective_eval_error(
+    hs071_initial_guess_fixture,
+    hs071_problem_instance_fixture,
+    hs071_definition_instance_fixture,
+    hs071_variable_lower_bounds_fixture,
+    hs071_variable_upper_bounds_fixture,
+    hs071_constraint_lower_bounds_fixture,
+    hs071_constraint_upper_bounds_fixture,
+):
+    class ObjectiveWithError:
+        def __init__(self):
+            self.n_eval_error = 0
+
+        def __call__(self, x):
+            if x[0] > 1.1:
+                self.n_eval_error += 1
+                raise cyipopt.CyIpoptEvaluationError()
+            return x[0] * x[3] * np.sum(x[0:3]) + x[2]
+
+    objective_with_error = ObjectiveWithError()
+
+    x0 = hs071_initial_guess_fixture
+    definition = hs071_definition_instance_fixture
+    definition.objective = objective_with_error
+    definition.intermediate = None
+
+    lb = hs071_variable_lower_bounds_fixture
+    ub = hs071_variable_upper_bounds_fixture
+    cl = hs071_constraint_lower_bounds_fixture
+    cu = hs071_constraint_upper_bounds_fixture
+
+    problem = _make_problem(definition, lb, ub, cl, cu)
+    # Note that the behavior tested here (success or failure of the solve when
+    # an evaluation error is raised in each callback) is documented in the
+    # CyIpoptEvaluationError class. If this behavior changes (e.g. Ipopt starts
+    # handling evaluation errors in the Jacobian), these tests will start to
+    # fail. We will need to (a) update these tests and (b) update the
+    # CyIpoptEvaluationError documentation, possibly with Ipopt version-specific
+    # behavior.
+    _solve_and_assert_correct(problem, x0)
+
+    assert objective_with_error.n_eval_error > 0
+
+
+def test_hs071_grad_obj_eval_error(
+    hs071_initial_guess_fixture,
+    hs071_problem_instance_fixture,
+    hs071_definition_instance_fixture,
+    hs071_variable_lower_bounds_fixture,
+    hs071_variable_upper_bounds_fixture,
+    hs071_constraint_lower_bounds_fixture,
+    hs071_constraint_upper_bounds_fixture,
+):
+    class GradObjWithError:
+        def __init__(self):
+            self.n_eval_error = 0
+
+        def __call__(self, x):
+            if x[0] > 1.1:
+                self.n_eval_error += 1
+                raise cyipopt.CyIpoptEvaluationError()
+            return np.array([
+                x[0] * x[3] + x[3] * np.sum(x[0:3]),
+                x[0] * x[3],
+                x[0] * x[3] + 1.0,
+                x[0] * np.sum(x[0:3]),
+            ])
+
+    gradient_with_error = GradObjWithError()
+
+    x0 = hs071_initial_guess_fixture
+    definition = hs071_definition_instance_fixture
+    definition.gradient = gradient_with_error
+    definition.intermediate = None
+
+    lb = hs071_variable_lower_bounds_fixture
+    ub = hs071_variable_upper_bounds_fixture
+    cl = hs071_constraint_lower_bounds_fixture
+    cu = hs071_constraint_upper_bounds_fixture
+
+    problem = _make_problem(definition, lb, ub, cl, cu)
+    # Solve fails when the evaluation error occurs in objective
+    # gradient evaluation.
+    _assert_solve_fails(problem, x0)
+
+    # Since we fail at the first evaluation error, we know we only encountered one.
+    assert gradient_with_error.n_eval_error == 1
+
+
+def test_hs071_constraints_eval_error(
+    hs071_initial_guess_fixture,
+    hs071_problem_instance_fixture,
+    hs071_definition_instance_fixture,
+    hs071_variable_lower_bounds_fixture,
+    hs071_variable_upper_bounds_fixture,
+    hs071_constraint_lower_bounds_fixture,
+    hs071_constraint_upper_bounds_fixture,
+):
+    class ConstraintsWithError:
+        def __init__(self):
+            self.n_eval_error = 0
+
+        def __call__(self, x):
+            if x[0] > 1.1:
+                self.n_eval_error += 1
+                raise cyipopt.CyIpoptEvaluationError()
+            return np.array((np.prod(x), np.dot(x, x)))
+
+    constraints_with_error = ConstraintsWithError()
+
+    x0 = hs071_initial_guess_fixture
+    definition = hs071_definition_instance_fixture
+    definition.constraints = constraints_with_error
+    definition.intermediate = None
+
+    lb = hs071_variable_lower_bounds_fixture
+    ub = hs071_variable_upper_bounds_fixture
+    cl = hs071_constraint_lower_bounds_fixture
+    cu = hs071_constraint_upper_bounds_fixture
+
+    problem = _make_problem(definition, lb, ub, cl, cu)
+    _solve_and_assert_correct(problem, x0)
+
+    assert constraints_with_error.n_eval_error > 0
+
+
+def test_hs071_jacobian_eval_error(
+    hs071_initial_guess_fixture,
+    hs071_problem_instance_fixture,
+    hs071_definition_instance_fixture,
+    hs071_variable_lower_bounds_fixture,
+    hs071_variable_upper_bounds_fixture,
+    hs071_constraint_lower_bounds_fixture,
+    hs071_constraint_upper_bounds_fixture,
+):
+    class JacobianWithError:
+        def __init__(self):
+            self.n_eval_error = 0
+
+        def __call__(self, x):
+            if x[0] > 1.1:
+                self.n_eval_error += 1
+                raise cyipopt.CyIpoptEvaluationError()
+            return np.concatenate((np.prod(x) / x, 2 * x))
+
+    jacobian_with_error = JacobianWithError()
+
+    x0 = hs071_initial_guess_fixture
+    definition = hs071_definition_instance_fixture
+    definition.jacobian = jacobian_with_error
+    definition.intermediate = None
+
+    lb = hs071_variable_lower_bounds_fixture
+    ub = hs071_variable_upper_bounds_fixture
+    cl = hs071_constraint_lower_bounds_fixture
+    cu = hs071_constraint_upper_bounds_fixture
+
+    problem = _make_problem(definition, lb, ub, cl, cu)
+    # Solve fails when the evaluation error occurs in constraint
+    # Jacobian evaluation.
+    _assert_solve_fails(problem, x0)
+
+    assert jacobian_with_error.n_eval_error == 1
+
+
+def test_hs071_hessian_eval_error(
+    hs071_initial_guess_fixture,
+    hs071_problem_instance_fixture,
+    hs071_definition_instance_fixture,
+    hs071_variable_lower_bounds_fixture,
+    hs071_variable_upper_bounds_fixture,
+    hs071_constraint_lower_bounds_fixture,
+    hs071_constraint_upper_bounds_fixture,
+):
+    class HessianWithError:
+        def __init__(self):
+            self.n_eval_error = 0
+
+        def __call__(self, x, lagrange, obj_factor):
+            if x[0] > 1.1:
+                self.n_eval_error += 1
+                raise cyipopt.CyIpoptEvaluationError()
+
+            H = obj_factor * np.array([
+                (2 * x[3], 0, 0, 0),
+                (x[3], 0, 0, 0),
+                (x[3], 0, 0, 0),
+                (2 * x[0] + x[1] + x[2], x[0], x[0], 0),
+            ])
+            H += lagrange[0] * np.array([
+                (0, 0, 0, 0),
+                (x[2] * x[3], 0, 0, 0),
+                (x[1] * x[3], x[0] * x[3], 0, 0),
+                (x[1] * x[2], x[0] * x[2], x[0] * x[1], 0),
+            ])
+            H += lagrange[1] * 2 * np.eye(4)
+            row, col = np.nonzero(np.tril(np.ones((4, 4))))
+            return H[row, col]
+
+    hessian_with_error = HessianWithError()
+
+    x0 = hs071_initial_guess_fixture
+    definition = hs071_definition_instance_fixture
+    definition.hessian = hessian_with_error
+    definition.intermediate = None
+
+    lb = hs071_variable_lower_bounds_fixture
+    ub = hs071_variable_upper_bounds_fixture
+    cl = hs071_constraint_lower_bounds_fixture
+    cu = hs071_constraint_upper_bounds_fixture
+
+    problem = _make_problem(definition, lb, ub, cl, cu)
+    # Solve fails when the evaluation error occurs in Lagrangian
+    # Hessian evaluation.
+    _assert_solve_fails(problem, x0)
+
+    assert hessian_with_error.n_eval_error == 1

docs/source/reference.rst

@@ -17,3 +17,5 @@ specifications may not be enough to give full guidelines on their uses.
 .. autofunction:: cyipopt.set_logging_level
 
 .. autofunction:: cyipopt.setLoggingLevel
+
+.. autoclass:: cyipopt.CyIpoptEvaluationError