Change gradient, jacobian, constraint_forward to return void

Values are now accessed directly from the problem struct:
- prob->constraint_values for constraint forward results
- prob->gradient_values for objective gradient
- prob->stacked_jac for constraint jacobian

Also updated Python bindings and tests to use the new API.

🤖 Generated with [Claude Code](https://claude.com/claude-code)

Co-Authored-By: Claude Opus 4.5 <noreply@anthropic.com>

Author: Transurgeon

include/problem.h

@@ -12,9 +12,11 @@ typedef struct problem
     int n_vars;
     int total_constraint_size;
 
-    /* Allocated by problem_init_derivatives */
+    /* Allocated by new_problem */
     double *constraint_values;
     double *gradient_values;
+
+    /* Allocated by problem_init_derivatives */
     CSR_Matrix *stacked_jac;
 } problem;
 
@@ -24,8 +26,8 @@ void problem_init_derivatives(problem *prob);
 void free_problem(problem *prob);
 
 double problem_objective_forward(problem *prob, const double *u);
-double *problem_constraint_forward(problem *prob, const double *u);
-double *problem_gradient(problem *prob, const double *u);
-CSR_Matrix *problem_jacobian(problem *prob, const double *u);
+void problem_constraint_forward(problem *prob, const double *u);
+void problem_gradient(problem *prob);
+void problem_jacobian(problem *prob);
 
 #endif

python/convert.py

@@ -151,11 +151,11 @@ def constraint_forward(self, u: np.ndarray) -> np.ndarray:
         """Evaluate constraints only. Returns constraint_values array."""
         return diffengine.problem_constraint_forward(self._capsule, u)
 
-    def gradient(self, u: np.ndarray) -> np.ndarray:
-        """Compute gradient of objective. Returns gradient array."""
-        return diffengine.problem_gradient(self._capsule, u)
+    def gradient(self) -> np.ndarray:
+        """Compute gradient of objective. Call objective_forward first. Returns gradient array."""
+        return diffengine.problem_gradient(self._capsule)
 
-    def jacobian(self, u: np.ndarray) -> sparse.csr_matrix:
-        """Compute jacobian of constraints. Returns scipy CSR matrix."""
-        data, indices, indptr, shape = diffengine.problem_jacobian(self._capsule, u)
+    def jacobian(self) -> sparse.csr_matrix:
+        """Compute jacobian of constraints. Call constraint_forward first. Returns scipy CSR matrix."""
+        data, indices, indptr, shape = diffengine.problem_jacobian(self._capsule)
         return sparse.csr_matrix((data, indices, indptr), shape=shape)

python/problem/constraint_forward.h

@@ -27,8 +27,8 @@ static PyObject *py_problem_constraint_forward(PyObject *self, PyObject *args)
         return NULL;
     }
 
-    double *constraint_vals =
-        problem_constraint_forward(prob, (const double *) PyArray_DATA(u_array));
+    problem_constraint_forward(prob, (const double *) PyArray_DATA(u_array));
+    Py_DECREF(u_array);
 
     PyObject *out = NULL;
     if (prob->total_constraint_size > 0)
@@ -37,10 +37,9 @@ static PyObject *py_problem_constraint_forward(PyObject *self, PyObject *args)
         out = PyArray_SimpleNew(1, &size, NPY_DOUBLE);
         if (!out)
         {
-            Py_DECREF(u_array);
             return NULL;
         }
-        memcpy(PyArray_DATA((PyArrayObject *) out), constraint_vals,
+        memcpy(PyArray_DATA((PyArrayObject *) out), prob->constraint_values,
                size * sizeof(double));
     }
     else
@@ -49,7 +48,6 @@ static PyObject *py_problem_constraint_forward(PyObject *self, PyObject *args)
         out = PyArray_SimpleNew(1, &size, NPY_DOUBLE);
     }
 
-    Py_DECREF(u_array);
     return out;
 }
 

python/problem/gradient.h

@@ -6,8 +6,7 @@
 static PyObject *py_problem_gradient(PyObject *self, PyObject *args)
 {
     PyObject *prob_capsule;
-    PyObject *u_obj;
-    if (!PyArg_ParseTuple(args, "OO", &prob_capsule, &u_obj))
+    if (!PyArg_ParseTuple(args, "O", &prob_capsule))
     {
         return NULL;
     }
@@ -20,25 +19,16 @@ static PyObject *py_problem_gradient(PyObject *self, PyObject *args)
         return NULL;
     }
 
-    PyArrayObject *u_array =
-        (PyArrayObject *) PyArray_FROM_OTF(u_obj, NPY_DOUBLE, NPY_ARRAY_IN_ARRAY);
-    if (!u_array)
-    {
-        return NULL;
-    }
-
-    double *grad = problem_gradient(prob, (const double *) PyArray_DATA(u_array));
+    problem_gradient(prob);
 
     npy_intp size = prob->n_vars;
     PyObject *out = PyArray_SimpleNew(1, &size, NPY_DOUBLE);
     if (!out)
     {
-        Py_DECREF(u_array);
         return NULL;
     }
-    memcpy(PyArray_DATA((PyArrayObject *) out), grad, size * sizeof(double));
+    memcpy(PyArray_DATA((PyArrayObject *) out), prob->gradient_values, size * sizeof(double));
 
-    Py_DECREF(u_array);
     return out;
 }
 

python/problem/jacobian.h

@@ -6,8 +6,7 @@
 static PyObject *py_problem_jacobian(PyObject *self, PyObject *args)
 {
     PyObject *prob_capsule;
-    PyObject *u_obj;
-    if (!PyArg_ParseTuple(args, "OO", &prob_capsule, &u_obj))
+    if (!PyArg_ParseTuple(args, "O", &prob_capsule))
     {
         return NULL;
     }
@@ -32,15 +31,9 @@ static PyObject *py_problem_jacobian(PyObject *self, PyObject *args)
         return Py_BuildValue("(OOO(ii))", data, indices, indptr, 0, prob->n_vars);
     }
 
-    PyArrayObject *u_array =
-        (PyArrayObject *) PyArray_FROM_OTF(u_obj, NPY_DOUBLE, NPY_ARRAY_IN_ARRAY);
-    if (!u_array)
-    {
-        return NULL;
-    }
-
-    CSR_Matrix *jac = problem_jacobian(prob, (const double *) PyArray_DATA(u_array));
+    problem_jacobian(prob);
 
+    CSR_Matrix *jac = prob->stacked_jac;
     npy_intp nnz = jac->nnz;
     npy_intp m_plus_1 = jac->m + 1;
 
@@ -53,15 +46,13 @@ static PyObject *py_problem_jacobian(PyObject *self, PyObject *args)
         Py_XDECREF(data);
         Py_XDECREF(indices);
         Py_XDECREF(indptr);
-        Py_DECREF(u_array);
         return NULL;
     }
 
     memcpy(PyArray_DATA((PyArrayObject *) data), jac->x, nnz * sizeof(double));
     memcpy(PyArray_DATA((PyArrayObject *) indices), jac->i, nnz * sizeof(int));
     memcpy(PyArray_DATA((PyArrayObject *) indptr), jac->p, m_plus_1 * sizeof(int));
 
-    Py_DECREF(u_array);
     return Py_BuildValue("(OOO(ii))", data, indices, indptr, jac->m, jac->n);
 }
 

python/tests/test_constrained.py

@@ -17,14 +17,14 @@ def test_single_constraint():
     prob = C_problem(problem)
 
     u = np.array([1.0, 2.0, 3.0])
-    prob.allocate(u)
+    prob.init_derivatives()
 
     # Test constraint_forward: constr.expr = -log(x)
     constraint_vals = prob.constraint_forward(u)
     assert np.allclose(constraint_vals, -np.log(u))
 
     # Test jacobian: d/dx(-log(x)) = -1/x
-    jac = prob.jacobian(u)
+    jac = prob.jacobian()
     expected_jac = np.diag(-1.0 / u)
     assert np.allclose(jac.toarray(), expected_jac)
 
@@ -42,15 +42,15 @@ def test_two_constraints():
     prob = C_problem(problem)
 
     u = np.array([1.0, 2.0])
-    prob.allocate(u)
+    prob.init_derivatives()
 
     # Test constraint_forward: [-log(u), -exp(u)]
     expected_constraint_vals = np.concatenate([-np.log(u), -np.exp(u)])
     constraint_vals = prob.constraint_forward(u)
     assert np.allclose(constraint_vals, expected_constraint_vals)
 
     # Test jacobian - stacked vertically
-    jac = prob.jacobian(u)
+    jac = prob.jacobian()
     assert jac.shape == (4, 2)
     expected_jac = np.vstack([np.diag(-1.0 / u), np.diag(-np.exp(u))])
     assert np.allclose(jac.toarray(), expected_jac)
@@ -70,7 +70,7 @@ def test_three_constraints_different_sizes():
     prob = C_problem(problem)
 
     u = np.array([1.0, 2.0, 3.0])
-    prob.allocate(u)
+    prob.init_derivatives()
 
     # Test constraint_forward
     expected_constraint_vals = np.concatenate([
@@ -82,7 +82,7 @@ def test_three_constraints_different_sizes():
     assert np.allclose(constraint_vals, expected_constraint_vals)
 
     # Test jacobian shape and values
-    jac = prob.jacobian(u)
+    jac = prob.jacobian()
     assert jac.shape == (7, 3)
     # First 3 rows: -diag(1/u), next 3 rows: -diag(exp(u)), last row: -1/u
     expected_jac = np.zeros((7, 3))
@@ -108,15 +108,15 @@ def test_multiple_variables():
     x_vals = np.array([1.0, 2.0])
     y_vals = np.array([0.5, 1.0])
     u = np.concatenate([x_vals, y_vals])
-    prob.allocate(u)
+    prob.init_derivatives()
 
     # Test constraint_forward
     expected_constraint_vals = np.concatenate([-np.log(x_vals), -np.exp(y_vals)])
     constraint_vals = prob.constraint_forward(u)
     assert np.allclose(constraint_vals, expected_constraint_vals)
 
     # Test jacobian
-    jac = prob.jacobian(u)
+    jac = prob.jacobian()
     assert jac.shape == (4, 4)
     expected_jac = np.zeros((4, 4))
     expected_jac[0, 0] = -1.0 / x_vals[0]
@@ -142,7 +142,7 @@ def test_larger_scale():
     prob = C_problem(problem)
 
     u = np.linspace(1.0, 5.0, n)
-    prob.allocate(u)
+    prob.init_derivatives()
 
     # Test constraint_forward
     expected_constraint_vals = np.concatenate([
@@ -155,7 +155,7 @@ def test_larger_scale():
     assert np.allclose(constraint_vals, expected_constraint_vals)
 
     # Test jacobian shape
-    jac = prob.jacobian(u)
+    jac = prob.jacobian()
     assert jac.shape == (n + n + 1 + 1, n)
 
 
@@ -169,16 +169,16 @@ def test_repeated_evaluations():
     prob = C_problem(problem)
 
     u1 = np.array([1.0, 2.0, 3.0])
-    prob.allocate(u1)
+    prob.init_derivatives()
 
     # First evaluation
     constraint_vals1 = prob.constraint_forward(u1)
-    jac1 = prob.jacobian(u1)
+    jac1 = prob.jacobian()
 
     # Second evaluation at different point
     u2 = np.array([2.0, 3.0, 4.0])
     constraint_vals2 = prob.constraint_forward(u2)
-    jac2 = prob.jacobian(u2)
+    jac2 = prob.jacobian()
 
     assert np.allclose(constraint_vals1, -np.exp(u1))
     assert np.allclose(constraint_vals2, -np.exp(u2))

python/tests/test_problem_native.py

@@ -58,7 +58,7 @@ def test_problem_gradient():
     diffengine.problem_init_derivatives(prob)
 
     diffengine.problem_objective_forward(prob, u)
-    grad = diffengine.problem_gradient(prob, u)
+    grad = diffengine.problem_gradient(prob)
     expected_grad = 1.0 / u
     assert np.allclose(grad, expected_grad)
 
@@ -76,7 +76,7 @@ def test_problem_jacobian():
     diffengine.problem_init_derivatives(prob)
 
     diffengine.problem_constraint_forward(prob, u)
-    data, indices, indptr, shape = diffengine.problem_jacobian(prob, u)
+    data, indices, indptr, shape = diffengine.problem_jacobian(prob)
     jac = sparse.csr_matrix((data, indices, indptr), shape=shape)
 
     expected_jac = np.diag(1.0 / u)
@@ -100,11 +100,11 @@ def test_problem_no_constraints():
     assert len(constraint_vals) == 0
 
     diffengine.problem_objective_forward(prob, u)
-    grad = diffengine.problem_gradient(prob, u)
+    grad = diffengine.problem_gradient(prob)
     assert np.allclose(grad, 1.0 / u)
 
     diffengine.problem_constraint_forward(prob, u)
-    data, indices, indptr, shape = diffengine.problem_jacobian(prob, u)
+    data, indices, indptr, shape = diffengine.problem_jacobian(prob)
     jac = sparse.csr_matrix((data, indices, indptr), shape=shape)
     assert jac.shape == (0, 3)
 
@@ -138,13 +138,13 @@ def test_problem_multiple_constraints():
 
     # Test gradient
     diffengine.problem_objective_forward(prob, u)
-    grad = diffengine.problem_gradient(prob, u)
+    grad = diffengine.problem_gradient(prob)
     expected_grad = 1.0 / u
     assert np.allclose(grad, expected_grad)
 
     # Test Jacobian
     diffengine.problem_constraint_forward(prob, u)
-    data, indices, indptr, shape = diffengine.problem_jacobian(prob, u)
+    data, indices, indptr, shape = diffengine.problem_jacobian(prob)
     jac = sparse.csr_matrix((data, indices, indptr), shape=shape)
 
     # Expected Jacobian: