Add get_current_iterate and get_current_violations methods to Problem class

cyipopt/cython/ipopt.pxd

@@ -167,3 +167,29 @@ cdef extern from "IpStdCInterface.h":
                     Number* mult_x_U,
                     UserDataPtr user_data
                     )
+
+    Bool GetIpoptCurrentIterate(
+                    IpoptProblem ipopt_problem,
+                    Bool scaled,
+                    Index n,
+                    Number* x,
+                    Number* z_L,
+                    Number* z_U,
+                    Index m,
+                    Number* g,
+                    Number* lambd
+                    )
+
+    Bool GetIpoptCurrentViolations(
+                    IpoptProblem ipopt_problem,
+                    Bool scaled,
+                    Index n,
+                    Number* x_L_violation,
+                    Number* x_U_violation,
+                    Number* compl_x_L,
+                    Number* compl_x_U,
+                    Number* grad_lag_x,
+                    Index m,
+                    Number* nlp_constraint_violation,
+                    Number* compl_g
+                    )

cyipopt/cython/ipopt_wrapper.pyx

@@ -654,6 +654,102 @@ cdef class Problem:
 
         return np_x, info
 
+    def get_current_iterate(self, scaled=False):
+        # Allocate arrays to hold the current iterate
+        cdef np.ndarray[DTYPEd_t, ndim=1] np_x
+        cdef np.ndarray[DTYPEd_t, ndim=1] np_mult_x_L
+        cdef np.ndarray[DTYPEd_t, ndim=1] np_mult_x_U
+        cdef np.ndarray[DTYPEd_t, ndim=1] np_g
+        cdef np.ndarray[DTYPEd_t, ndim=1] np_mult_g
+        np_x = np.zeros((self.__n,), dtype=DTYPEd).flatten()
+        np_mult_x_L = np.zeros((self.__n,), dtype=DTYPEd).flatten()
+        np_mult_x_U = np.zeros((self.__n,), dtype=DTYPEd).flatten()
+        np_g = np.zeros((self.__m,), dtype=DTYPEd).flatten()
+        np_mult_g = np.zeros((self.__m,), dtype=DTYPEd).flatten()
+
+        # Cast to C data types
+        x = <Number*>np_x.data
+        mult_x_L = <Number*>np_mult_x_L.data
+        mult_x_U = <Number*>np_mult_x_U.data
+        g = <Number*>np_g.data
+        mult_g = <Number*>np_mult_g.data
+
+        # NOTE: GetIpoptCurrentIterate can *only* be called during an
+        # intermediate callback (otherwise __nlp->tnlp->ip_data_ is NULL)
+        # TODO: Either catch error or avoid calling if we are not in an
+        # intermediate callback
+        ret = GetIpoptCurrentIterate(
+            self.__nlp,
+            scaled,
+            self.__n,
+            x,
+            mult_x_L,
+            mult_x_U,
+            self.__m,
+            g,
+            mult_g,
+        )
+
+        # Return values to user
+        # - Is another data type (e.g. dict, namedtuple) more appropriate than
+        #   simply a tuple?
+        # - Should `ret` be returned here?
+        return (np_x, np_mult_x_L, np_mult_x_U, np_g, np_mult_g)
+
+    def get_current_violations(self, scaled=False):
+        # Allocate arrays to hold current violations
+        cdef np.ndarray[DTYPEd_t, ndim=1] np_x_L_viol
+        cdef np.ndarray[DTYPEd_t, ndim=1] np_x_U_viol
+        cdef np.ndarray[DTYPEd_t, ndim=1] np_compl_x_L
+        cdef np.ndarray[DTYPEd_t, ndim=1] np_compl_x_U
+        cdef np.ndarray[DTYPEd_t, ndim=1] np_grad_lag_x
+        cdef np.ndarray[DTYPEd_t, ndim=1] np_g_viol
+        cdef np.ndarray[DTYPEd_t, ndim=1] np_compl_g
+        np_x_L_viol = np.zeros((self.__n,), dtype=DTYPEd).flatten()
+        np_x_U_viol = np.zeros((self.__n,), dtype=DTYPEd).flatten()
+        np_compl_x_L = np.zeros((self.__n,), dtype=DTYPEd).flatten()
+        np_compl_x_U = np.zeros((self.__n,), dtype=DTYPEd).flatten()
+        np_grad_lag_x = np.zeros((self.__n,), dtype=DTYPEd).flatten()
+        np_g_viol = np.zeros((self.__m,), dtype=DTYPEd).flatten()
+        np_compl_g = np.zeros((self.__m,), dtype=DTYPEd).flatten()
+
+        # Cast to C data types
+        x_L_viol = <Number*>np_x_L_viol.data
+        x_U_viol = <Number*>np_x_U_viol.data
+        compl_x_L = <Number*>np_compl_x_L.data
+        compl_x_U = <Number*>np_compl_x_U.data
+        grad_lag_x = <Number*>np_grad_lag_x.data
+        g_viol = <Number*>np_g_viol.data
+        compl_g = <Number*>np_compl_g.data
+
+        # NOTE: GetIpoptCurrentViolations can *only* be called during an
+        # intermediate callback (otherwise __nlp->tnlp->ip_data_ is NULL)
+        # TODO: Either catch error or avoid calling if we are not in an
+        # intermediate callback
+        ret = GetIpoptCurrentViolations(
+            self.__nlp,
+            scaled,
+            self.__n,
+            x_L_viol,
+            x_U_viol,
+            compl_x_L,
+            compl_x_U,
+            grad_lag_x,
+            self.__m,
+            g_viol,
+            compl_g,
+        )
+
+        return (
+            np_x_L_viol,
+            np_x_U_viol,
+            np_compl_x_L,
+            np_compl_x_U,
+            np_grad_lag_x,
+            np_g_viol,
+            np_compl_g,
+        )
+
 
 #
 # Callback functions

cyipopt/tests/conftest.py

@@ -111,7 +111,7 @@ def intermediate(*args):
 
 
 @pytest.fixture()
-def hs071_defintion_instance_fixture(hs071_objective_fixture,
+def hs071_definition_instance_fixture(hs071_objective_fixture,
                                      hs071_gradient_fixture,
                                      hs071_constraints_fixture,
                                      hs071_jacobian_fixture,
@@ -175,15 +175,15 @@ def hs071_constraint_upper_bounds_fixture():
 
 
 @pytest.fixture()
-def hs071_problem_instance_fixture(hs071_defintion_instance_fixture,
+def hs071_problem_instance_fixture(hs071_definition_instance_fixture,
                                    hs071_initial_guess_fixture,
                                    hs071_variable_lower_bounds_fixture,
                                    hs071_variable_upper_bounds_fixture,
                                    hs071_constraint_lower_bounds_fixture,
                                    hs071_constraint_upper_bounds_fixture,
                                    ):
     """Return a default cyipopt.Problem instance of the hs071 test problem."""
-    problem_definition = hs071_defintion_instance_fixture
+    problem_definition = hs071_definition_instance_fixture
     x0 = hs071_initial_guess_fixture
     lb = hs071_variable_lower_bounds_fixture
     ub = hs071_variable_upper_bounds_fixture

cyipopt/tests/unit/test_deprecations.py

@@ -22,7 +22,7 @@ def test_ipopt_import_deprecation():
         import ipopt
 
 
-def test_non_pep8_class_name_deprecation(hs071_defintion_instance_fixture,
+def test_non_pep8_class_name_deprecation(hs071_definition_instance_fixture,
                                          hs071_initial_guess_fixture,
                                          hs071_variable_lower_bounds_fixture,
                                          hs071_variable_upper_bounds_fixture,
@@ -36,7 +36,7 @@ def test_non_pep8_class_name_deprecation(hs071_defintion_instance_fixture,
     with pytest.warns(FutureWarning, match=expected_warning_msg):
         _ = cyipopt.problem(n=len(hs071_initial_guess_fixture),
                             m=len(hs071_constraint_lower_bounds_fixture),
-                            problem_obj=hs071_defintion_instance_fixture,
+                            problem_obj=hs071_definition_instance_fixture,
                             lb=hs071_variable_lower_bounds_fixture,
                             ub=hs071_variable_upper_bounds_fixture,
                             cl=hs071_constraint_lower_bounds_fixture,

cyipopt/tests/unit/test_ipopt_funcs.py

@@ -0,0 +1,213 @@
+import numpy as np
+import pytest
+
+import cyipopt
+
+@pytest.mark.skipif(True, reason="This segfaults. Ideally, it fails gracefully")
+def test_get_iterate_uninit(hs071_problem_instance_fixture):
+    """Test that we can call get_current_iterate on an uninitialized problem
+    """
+    nlp = hs071_problem_instance_fixture
+    x, zL, zU, g, lam = nlp.get_current_iterate()
+
+
+@pytest.mark.skipif(True, reason="This also segfaults")
+def test_get_iterate_postsolve(
+    hs071_initial_guess_fixture,
+    hs071_problem_instance_fixture,
+):
+    x0 = hs071_initial_guess_fixture
+    nlp = hs071_problem_instance_fixture
+    x, info = nlp.solve(x0)
+
+    x, zL, zU, g, lam = nlp.get_current_iterate()
+    expected_x = np.array([1.0, 4.74299964, 3.82114998, 1.37940829])
+    np.testing.assert_allclose(x, expected_x)
+
+
+@pytest.mark.skipif(True, reason="Segfaults")
+def test_get_violations_uninit(hs071_problem_instance_fixture):
+    nlp = hs071_problem_instance_fixture
+    x, zL, zU, g, lam = nlp.get_current_violations()
+
+
+@pytest.mark.skipif(True, reason="Segfaults")
+def test_get_violations_postsolve(
+    hs071_initial_guess_fixture,
+    hs071_problem_instance_fixture,
+):
+    x0 = hs071_initial_guess_fixture
+    nlp = hs071_problem_instance_fixture
+    x, info = nlp.solve(x0)
+
+    x, zL, zU, g, lam = nlp.get_current_violations()
+    expected_x = np.array([1.0, 4.74299964, 3.82114998, 1.37940829])
+    np.testing.assert_allclose(x, expected_x)
+
+
+def test_get_iterate_hs071(
+    hs071_initial_guess_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,
+):
+    x0 = hs071_initial_guess_fixture
+    lb = hs071_variable_lower_bounds_fixture
+    ub = hs071_variable_upper_bounds_fixture
+    cl = hs071_constraint_lower_bounds_fixture
+    cu = hs071_constraint_upper_bounds_fixture
+    n = len(x0)
+    m = len(cl)
+
+    problem_definition = hs071_definition_instance_fixture
+
+    #
+    # Define a callback that uses some "global" information to call
+    # get_current_iterate and store the result
+    #
+    x_iterates = []
+    def intermediate(
+        alg_mod,
+        iter_count,
+        obj_value,
+        inf_pr,
+        inf_du,
+        mu,
+        d_norm,
+        regularization_size,
+        alpha_du,
+        alpha_pr,
+        ls_trials,
+    ):
+        # CyIpopt's C wapper expects a callback with this signature. If we
+        # implemented this as a method on problem_definition, we could store
+        # and access global information on self.
+
+        # This callback must be defined before constructing the Problem, but can
+        # be defined after (or as part of) problem_definition. If we attach the
+        # Problem to the "definition", then we can call get_current_iterate
+        # from this callback.
+        iterate = problem_definition.nlp.get_current_iterate(scaled=False)
+        x, zL, zU, g, lam = iterate
+        x_iterates.append(x)
+
+        # Hack so we may get the number of iterations after the solve
+        problem_definition.iter_count = iter_count
+
+    # Replace "intermediate" attribute with our callback, which knows
+    # about the "Problem", and therefore can call get_current_iterate.
+    problem_definition.intermediate = intermediate
+
+    nlp = cyipopt.Problem(
+        n=n,
+        m=m,
+        problem_obj=problem_definition,
+        lb=lb,
+        ub=ub,
+        cl=cl,
+        cu=cu,
+    )
+    # Add nlp (the "Problem") as an attribute on our "problem definition".
+    # This way we can call methods on the Problem, like get_current_iterate,
+    # during the solve.
+    problem_definition.nlp = nlp
+
+    # Disable bound push to make testing easier
+    nlp.add_option("bound_push", 1e-9)
+    x, info = nlp.solve(x0)
+
+    # Assert correct solution
+    expected_x = np.array([1.0, 4.74299964, 3.82114998, 1.37940829])
+    np.testing.assert_allclose(x, expected_x)
+
+    #
+    # Assert some very basic information about the collected primal iterates
+    #
+    assert len(x_iterates) == (1 + problem_definition.iter_count)
+
+    # These could be different due to bound_push (and scaling)
+    np.testing.assert_allclose(x_iterates[0], x0)
+
+    # These could be different due to honor_original_bounds (and scaling)
+    np.testing.assert_allclose(x_iterates[-1], x)
+
+
+def test_get_violations_hs071(
+    hs071_initial_guess_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,
+):
+    x0 = hs071_initial_guess_fixture
+    lb = hs071_variable_lower_bounds_fixture
+    ub = hs071_variable_upper_bounds_fixture
+    cl = hs071_constraint_lower_bounds_fixture
+    cu = hs071_constraint_upper_bounds_fixture
+    n = len(x0)
+    m = len(cl)
+
+    problem_definition = hs071_definition_instance_fixture
+
+    pr_violations = []
+    du_violations = []
+    def intermediate(
+        alg_mod,
+        iter_count,
+        obj_value,
+        inf_pr,
+        inf_du,
+        mu,
+        d_norm,
+        regularization_size,
+        alpha_du,
+        alpha_pr,
+        ls_trials,
+    ):
+        violations = problem_definition.nlp.get_current_violations(scaled=True)
+        (
+            xL_viol, xU_viol, xL_compl, xU_compl, grad_lag, g_viol, g_compl
+        ) = violations
+        pr_violations.append(g_viol)
+        du_violations.append(grad_lag)
+
+        # Hack so we may get the number of iterations after the solve
+        problem_definition.iter_count = iter_count
+
+    problem_definition.intermediate = intermediate
+    nlp = cyipopt.Problem(
+        n=n,
+        m=m,
+        problem_obj=problem_definition,
+        lb=lb,
+        ub=ub,
+        cl=cl,
+        cu=cu,
+    )
+    problem_definition.nlp = nlp
+
+    nlp.add_option("tol", 1e-8)
+    x, info = nlp.solve(x0)
+
+    # Assert correct solution
+    expected_x = np.array([1.0, 4.74299964, 3.82114998, 1.37940829])
+    np.testing.assert_allclose(x, expected_x)
+
+    #
+    # Assert some very basic information about the collected violations
+    #
+    assert len(pr_violations) == (1 + problem_definition.iter_count)
+    assert len(du_violations) == (1 + problem_definition.iter_count)
+
+    #
+    # With atol=1e-8, this check fails. This differs from what I see in the
+    # Ipopt log, where inf_pr is 1.77e-11 at the final iteration. I see
+    # final primal violations: [2.455637e-07, 1.770672e-11]
+    # Not sure if a bug or not...
+    #
+    np.testing.assert_allclose(pr_violations[-1], np.zeros(m), atol=1e-6)
+
+    np.testing.assert_allclose(du_violations[-1], np.zeros(n), atol=1e-8)