Adds convert function from cvxpy problem to C

docs/problem_struct_design.md

@@ -0,0 +1,205 @@
+# Design: C Problem Struct
+
+## Summary
+
+Create a native C `problem` struct that encapsulates objective + constraints, with methods for:
+- `forward(u)` - evaluate objective and all constraints
+- `gradient(u)` - return objective gradient (jacobian.T)
+- `jacobian(u)` - return single stacked CSR matrix of constraint jacobians
+
+## Files to Create/Modify
+
+1. **`include/problem.h`** - New header defining problem struct
+2. **`src/problem.c`** - Implementation
+3. **`python/bindings.c`** - Python bindings for problem
+4. **`python/convert.py`** - Update to return problem capsule
+5. **`tests/problem/test_problem.h`** - C tests
+
+---
+
+## Step 1: Create `include/problem.h`
+
+```c
+#ifndef PROBLEM_H
+#define PROBLEM_H
+
+#include "expr.h"
+#include "utils/CSR_Matrix.h"
+
+typedef struct problem
+{
+    expr *objective;           /* Objective expression (scalar) */
+    expr **constraints;        /* Array of constraint expressions */
+    int n_constraints;
+    int n_vars;
+    int total_constraint_size; /* Sum of all constraint sizes */
+
+    /* Pre-allocated storage */
+    double *constraint_values;
+    double *gradient_values;   /* Dense gradient array */
+    CSR_Matrix *stacked_jac;
+} problem;
+
+problem *new_problem(expr *objective, expr **constraints, int n_constraints);
+void free_problem(problem *prob);
+double problem_forward(problem *prob, const double *u);
+double *problem_gradient(problem *prob, const double *u);
+CSR_Matrix *problem_jacobian(problem *prob, const double *u);
+
+#endif
+```
+
+---
+
+## Step 2: Create `src/problem.c`
+
+Key functions:
+
+### `new_problem`
+- Retain (increment refcount) on objective and all constraints
+- Compute `total_constraint_size = sum(constraints[i]->size)`
+- Pre-allocate `constraint_values` and `gradient_values` arrays
+
+### `free_problem`
+- Call `free_expr` on objective and all constraints (decrements refcount)
+- Free allocated arrays and stacked_jac
+
+### `problem_forward`
+```c
+double problem_forward(problem *prob, const double *u)
+{
+    prob->objective->forward(prob->objective, u);
+    double obj_val = prob->objective->value[0];
+
+    int offset = 0;
+    for (int i = 0; i < prob->n_constraints; i++)
+    {
+        expr *c = prob->constraints[i];
+        c->forward(c, u);
+        memcpy(prob->constraint_values + offset, c->value, c->size * sizeof(double));
+        offset += c->size;
+    }
+    return obj_val;
+}
+```
+
+### `problem_gradient`
+- Run forward pass on objective
+- Call jacobian_init + eval_jacobian
+- Objective jacobian is 1 x n_vars row vector
+- Copy sparse row to dense `gradient_values` array
+- Return pointer to internal array
+
+### `problem_jacobian`
+- Forward + jacobian for each constraint
+- Stack CSR matrices vertically:
+  - Total rows = `total_constraint_size`
+  - Copy row pointers with offset, copy col indices and values
+- Lazy allocate/reallocate `stacked_jac` based on total nnz
+
+---
+
+## Step 3: Update `python/bindings.c`
+
+Add capsule and functions:
+
+```c
+#define PROBLEM_CAPSULE_NAME "DNLP_PROBLEM"
+
+static void problem_capsule_destructor(PyObject *capsule) { ... }
+
+static PyObject *py_make_problem(PyObject *self, PyObject *args)
+{
+    PyObject *obj_capsule, *constraints_list;
+    // Parse objective capsule and list of constraint capsules
+    // Extract expr* pointers, call new_problem
+    // Return PyCapsule
+}
+
+static PyObject *py_problem_forward(PyObject *self, PyObject *args)
+{
+    // Returns: (obj_value, constraint_values_array)
+}
+
+static PyObject *py_problem_gradient(PyObject *self, PyObject *args)
+{
+    // Returns: numpy array of size n_vars
+}
+
+static PyObject *py_problem_jacobian(PyObject *self, PyObject *args)
+{
+    // Returns: (data, indices, indptr, (m, n)) for scipy CSR
+}
+```
+
+Add to DNLPMethods:
+```c
+{"make_problem", py_make_problem, METH_VARARGS, "Create problem"},
+{"problem_forward", py_problem_forward, METH_VARARGS, "Evaluate problem"},
+{"problem_gradient", py_problem_gradient, METH_VARARGS, "Compute gradient"},
+{"problem_jacobian", py_problem_jacobian, METH_VARARGS, "Compute constraint jacobian"},
+```
+
+---
+
+## Step 4: Update `python/convert.py`
+
+```python
+def convert_problem(problem: cp.Problem):
+    """Convert CVXPY Problem to C problem struct."""
+    var_dict = build_variable_dict(problem.variables())
+
+    c_objective = _convert_expr(problem.objective.expr, var_dict)
+    c_constraints = [_convert_expr(c.expr, var_dict) for c in problem.constraints]
+
+    return diffengine.make_problem(c_objective, c_constraints)
+
+
+class Problem:
+    """Wrapper for C problem struct with clean Python API."""
+
+    def __init__(self, cvxpy_problem: cp.Problem):
+        self._capsule = convert_problem(cvxpy_problem)
+
+    def forward(self, u: np.ndarray) -> tuple[float, np.ndarray]:
+        return diffengine.problem_forward(self._capsule, u)
+
+    def gradient(self, u: np.ndarray) -> np.ndarray:
+        return diffengine.problem_gradient(self._capsule, u)
+
+    def jacobian(self, u: np.ndarray) -> sparse.csr_matrix:
+        data, indices, indptr, shape = diffengine.problem_jacobian(self._capsule, u)
+        return sparse.csr_matrix((data, indices, indptr), shape=shape)
+```
+
+---
+
+## Step 5: Add Tests
+
+### C tests in `tests/problem/test_problem.h`:
+- `test_problem_forward` - verify objective and constraint values
+- `test_problem_gradient` - verify gradient matches manual calculation
+- `test_problem_jacobian_stacking` - verify stacked matrix structure
+
+### Python tests in `convert.py`:
+- `test_problem_forward` - compare with numpy
+- `test_problem_gradient` - gradient of sum(log(x)) = 1/x
+- `test_problem_jacobian` - verify stacked jacobian shape and values
+
+---
+
+## Implementation Order
+
+1. Create `include/problem.h`
+2. Create `src/problem.c` with new_problem, free_problem, problem_forward
+3. Add problem_gradient and problem_jacobian
+4. Add Python bindings to `bindings.c`
+5. Rebuild: `cmake --build build`
+6. Update `convert.py` with Problem class
+7. Add and run tests
+
+## Key Design Notes
+
+- **Memory**: Uses expr refcounting - new_problem retains, free_problem releases
+- **Efficiency**: Pre-allocates arrays, lazy-allocates stacked jacobian
+- **API**: Returns internal pointers (caller should NOT free)

python/bindings.c

@@ -10,9 +10,13 @@
 // Capsule name for expr* pointers
 #define EXPR_CAPSULE_NAME "DNLP_EXPR"
 
+static int numpy_initialized = 0;
+
 static int ensure_numpy(void)
 {
-    import_array();
+    if (numpy_initialized) return 0;
+    import_array1(-1);
+    numpy_initialized = 1;
     return 0;
 }
 
@@ -65,6 +69,77 @@ static PyObject *py_make_log(PyObject *self, PyObject *args)
     return PyCapsule_New(node, EXPR_CAPSULE_NAME, expr_capsule_destructor);
 }
 
+static PyObject *py_make_exp(PyObject *self, PyObject *args)
+{
+    PyObject *child_capsule;
+    if (!PyArg_ParseTuple(args, "O", &child_capsule))
+    {
+        return NULL;
+    }
+    expr *child = (expr *) PyCapsule_GetPointer(child_capsule, EXPR_CAPSULE_NAME);
+    if (!child)
+    {
+        PyErr_SetString(PyExc_ValueError, "invalid child capsule");
+        return NULL;
+    }
+
+    expr *node = new_exp(child);
+    if (!node)
+    {
+        PyErr_SetString(PyExc_RuntimeError, "failed to create exp node");
+        return NULL;
+    }
+    return PyCapsule_New(node, EXPR_CAPSULE_NAME, expr_capsule_destructor);
+}
+
+static PyObject *py_make_add(PyObject *self, PyObject *args)
+{
+    PyObject *left_capsule, *right_capsule;
+    if (!PyArg_ParseTuple(args, "OO", &left_capsule, &right_capsule))
+    {
+        return NULL;
+    }
+    expr *left = (expr *) PyCapsule_GetPointer(left_capsule, EXPR_CAPSULE_NAME);
+    expr *right = (expr *) PyCapsule_GetPointer(right_capsule, EXPR_CAPSULE_NAME);
+    if (!left || !right)
+    {
+        PyErr_SetString(PyExc_ValueError, "invalid child capsule");
+        return NULL;
+    }
+
+    expr *node = new_add(left, right);
+    if (!node)
+    {
+        PyErr_SetString(PyExc_RuntimeError, "failed to create add node");
+        return NULL;
+    }
+    return PyCapsule_New(node, EXPR_CAPSULE_NAME, expr_capsule_destructor);
+}
+
+static PyObject *py_make_sum(PyObject *self, PyObject *args)
+{
+    PyObject *child_capsule;
+    int axis;
+    if (!PyArg_ParseTuple(args, "Oi", &child_capsule, &axis))
+    {
+        return NULL;
+    }
+    expr *child = (expr *) PyCapsule_GetPointer(child_capsule, EXPR_CAPSULE_NAME);
+    if (!child)
+    {
+        PyErr_SetString(PyExc_ValueError, "invalid child capsule");
+        return NULL;
+    }
+
+    expr *node = new_sum(child, axis);
+    if (!node)
+    {
+        PyErr_SetString(PyExc_RuntimeError, "failed to create sum node");
+        return NULL;
+    }
+    return PyCapsule_New(node, EXPR_CAPSULE_NAME, expr_capsule_destructor);
+}
+
 static PyObject *py_forward(PyObject *self, PyObject *args)
 {
     PyObject *node_capsule;
@@ -103,10 +178,73 @@ static PyObject *py_forward(PyObject *self, PyObject *args)
     return out;
 }
 
+static PyObject *py_jacobian(PyObject *self, PyObject *args)
+{
+    PyObject *node_capsule;
+    PyObject *u_obj;
+    if (!PyArg_ParseTuple(args, "OO", &node_capsule, &u_obj))
+    {
+        return NULL;
+    }
+
+    expr *node = (expr *) PyCapsule_GetPointer(node_capsule, EXPR_CAPSULE_NAME);
+    if (!node)
+    {
+        PyErr_SetString(PyExc_ValueError, "invalid node capsule");
+        return NULL;
+    }
+
+    PyArrayObject *u_array =
+        (PyArrayObject *) PyArray_FROM_OTF(u_obj, NPY_DOUBLE, NPY_ARRAY_IN_ARRAY);
+    if (!u_array)
+    {
+        return NULL;
+    }
+
+    // Run forward pass first (required before jacobian)
+    node->forward(node, (const double *) PyArray_DATA(u_array));
+
+    // Initialize and evaluate jacobian
+    node->jacobian_init(node);
+    node->eval_jacobian(node);
+
+    CSR_Matrix *jac = node->jacobian;
+
+    // Create numpy arrays for CSR components
+    npy_intp nnz = jac->nnz;
+    npy_intp m_plus_1 = jac->m + 1;
+
+    PyObject *data = PyArray_SimpleNew(1, &nnz, NPY_DOUBLE);
+    PyObject *indices = PyArray_SimpleNew(1, &nnz, NPY_INT32);
+    PyObject *indptr = PyArray_SimpleNew(1, &m_plus_1, NPY_INT32);
+
+    if (!data || !indices || !indptr)
+    {
+        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 tuple: (data, indices, indptr, shape)
+    return Py_BuildValue("(OOO(ii))", data, indices, indptr, jac->m, jac->n);
+}
+
 static PyMethodDef DNLPMethods[] = {
     {"make_variable", py_make_variable, METH_VARARGS, "Create variable node"},
     {"make_log", py_make_log, METH_VARARGS, "Create log node"},
+    {"make_exp", py_make_exp, METH_VARARGS, "Create exp node"},
+    {"make_add", py_make_add, METH_VARARGS, "Create add node"},
+    {"make_sum", py_make_sum, METH_VARARGS, "Create sum node"},
     {"forward", py_forward, METH_VARARGS, "Run forward pass and return values"},
+    {"jacobian", py_jacobian, METH_VARARGS, "Compute jacobian and return CSR components"},
     {NULL, NULL, 0, NULL}};
 
 static struct PyModuleDef dnlp_module = {PyModuleDef_HEAD_INIT, "DNLP_diff_engine",