diff --git a/CMakeLists.txt b/CMakeLists.txt
index 5d938dc..24f9139 100644
--- a/CMakeLists.txt
+++ b/CMakeLists.txt
@@ -61,25 +61,6 @@ endif()
 # Enable position-independent code for shared library compatibility
 set_property(TARGET dnlp_diff PROPERTY POSITION_INDEPENDENT_CODE ON)
 
-# =============================================================================
-# Python bindings (built when using scikit-build-core)
-# =============================================================================
-#if(SKBUILD)
-#    find_package(Python3 REQUIRED COMPONENTS Interpreter Development.Module NumPy)
-#
-#    # Create Python extension module
-#    Python3_add_library(_core MODULE python/bindings.c)
-#    target_include_directories(_core PRIVATE
-#        ${PROJECT_SOURCE_DIR}/include
-#        ${PROJECT_SOURCE_DIR}/python
-#        ${Python3_NumPy_INCLUDE_DIRS}
-#    )
-#    target_link_libraries(_core PRIVATE dnlp_diff)
-#
-#    # Install to the package directory
-#    install(TARGETS _core LIBRARY DESTINATION dnlp_diff_engine)
-#endif()
-
 # =============================================================================
 # C tests (only for standalone builds)
 # =============================================================================
diff --git a/python/atoms/add.h b/python/atoms/add.h
deleted file mode 100644
index 2f62a7a..0000000
--- a/python/atoms/add.h
+++ /dev/null
@@ -1,31 +0,0 @@
-#ifndef ATOM_ADD_H
-#define ATOM_ADD_H
-
-#include "common.h"
-
-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;
-    }
-    expr_retain(node); /* Capsule owns a reference */
-    return PyCapsule_New(node, EXPR_CAPSULE_NAME, expr_capsule_destructor);
-}
-
-#endif /* ATOM_ADD_H */
diff --git a/python/atoms/asinh.h b/python/atoms/asinh.h
deleted file mode 100644
index a09934e..0000000
--- a/python/atoms/asinh.h
+++ /dev/null
@@ -1,31 +0,0 @@
-#ifndef ATOM_ASINH_H
-#define ATOM_ASINH_H
-
-#include "common.h"
-#include "elementwise_univariate.h"
-
-static PyObject *py_make_asinh(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_asinh(child);
-    if (!node)
-    {
-        PyErr_SetString(PyExc_RuntimeError, "failed to create asinh node");
-        return NULL;
-    }
-    expr_retain(node); /* Capsule owns a reference */
-    return PyCapsule_New(node, EXPR_CAPSULE_NAME, expr_capsule_destructor);
-}
-
-#endif /* ATOM_ASINH_H */
diff --git a/python/atoms/atanh.h b/python/atoms/atanh.h
deleted file mode 100644
index fd0568d..0000000
--- a/python/atoms/atanh.h
+++ /dev/null
@@ -1,31 +0,0 @@
-#ifndef ATOM_ATANH_H
-#define ATOM_ATANH_H
-
-#include "common.h"
-#include "elementwise_univariate.h"
-
-static PyObject *py_make_atanh(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_atanh(child);
-    if (!node)
-    {
-        PyErr_SetString(PyExc_RuntimeError, "failed to create atanh node");
-        return NULL;
-    }
-    expr_retain(node); /* Capsule owns a reference */
-    return PyCapsule_New(node, EXPR_CAPSULE_NAME, expr_capsule_destructor);
-}
-
-#endif /* ATOM_ATANH_H */
diff --git a/python/atoms/broadcast.h b/python/atoms/broadcast.h
deleted file mode 100644
index 10a5e1e..0000000
--- a/python/atoms/broadcast.h
+++ /dev/null
@@ -1,33 +0,0 @@
-#ifndef ATOM_BROADCAST_H
-#define ATOM_BROADCAST_H
-
-#include "common.h"
-
-static PyObject *py_make_broadcast(PyObject *self, PyObject *args)
-{
-    PyObject *child_capsule;
-    int d1, d2;
-
-    if (!PyArg_ParseTuple(args, "Oii", &child_capsule, &d1, &d2))
-    {
-        return NULL;
-    }
-
-    expr *child = (expr *) PyCapsule_GetPointer(child_capsule, EXPR_CAPSULE_NAME);
-    if (!child)
-    {
-        return NULL;
-    }
-
-    expr *node = new_broadcast(child, d1, d2);
-    if (!node)
-    {
-        PyErr_SetString(PyExc_RuntimeError, "failed to create broadcast node");
-        return NULL;
-    }
-
-    expr_retain(node);
-    return PyCapsule_New(node, EXPR_CAPSULE_NAME, expr_capsule_destructor);
-}
-
-#endif /* ATOM_BROADCAST_H */
diff --git a/python/atoms/common.h b/python/atoms/common.h
deleted file mode 100644
index 83101dd..0000000
--- a/python/atoms/common.h
+++ /dev/null
@@ -1,23 +0,0 @@
-#ifndef ATOMS_COMMON_H
-#define ATOMS_COMMON_H
-
-#define PY_SSIZE_T_CLEAN
-#include <Python.h>
-#include <numpy/arrayobject.h>
-
-#include "affine.h"
-#include "elementwise_univariate.h"
-#include "expr.h"
-
-#define EXPR_CAPSULE_NAME "DNLP_EXPR"
-
-static void expr_capsule_destructor(PyObject *capsule)
-{
-    expr *node = (expr *) PyCapsule_GetPointer(capsule, EXPR_CAPSULE_NAME);
-    if (node)
-    {
-        free_expr(node);
-    }
-}
-
-#endif /* ATOMS_COMMON_H */
diff --git a/python/atoms/const_scalar_mult.h b/python/atoms/const_scalar_mult.h
deleted file mode 100644
index 1a83e5a..0000000
--- a/python/atoms/const_scalar_mult.h
+++ /dev/null
@@ -1,36 +0,0 @@
-#ifndef ATOM_CONST_SCALAR_MULT_H
-#define ATOM_CONST_SCALAR_MULT_H
-
-#include "bivariate.h"
-#include "common.h"
-
-/* Constant scalar multiplication: a * f(x) where a is a constant double */
-static PyObject *py_make_const_scalar_mult(PyObject *self, PyObject *args)
-{
-    PyObject *child_capsule;
-    double a;
-
-    if (!PyArg_ParseTuple(args, "Od", &child_capsule, &a))
-    {
-        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_const_scalar_mult(a, child);
-    if (!node)
-    {
-        PyErr_SetString(PyExc_RuntimeError,
-                        "failed to create const_scalar_mult node");
-        return NULL;
-    }
-    expr_retain(node); /* Capsule owns a reference */
-    return PyCapsule_New(node, EXPR_CAPSULE_NAME, expr_capsule_destructor);
-}
-
-#endif /* ATOM_CONST_SCALAR_MULT_H */
diff --git a/python/atoms/const_vector_mult.h b/python/atoms/const_vector_mult.h
deleted file mode 100644
index aa1f83e..0000000
--- a/python/atoms/const_vector_mult.h
+++ /dev/null
@@ -1,59 +0,0 @@
-#ifndef ATOM_CONST_VECTOR_MULT_H
-#define ATOM_CONST_VECTOR_MULT_H
-
-#include "bivariate.h"
-#include "common.h"
-
-/* Constant vector elementwise multiplication: a ∘ f(x) where a is a constant vector
- */
-static PyObject *py_make_const_vector_mult(PyObject *self, PyObject *args)
-{
-    PyObject *child_capsule;
-    PyObject *a_obj;
-
-    if (!PyArg_ParseTuple(args, "OO", &child_capsule, &a_obj))
-    {
-        return NULL;
-    }
-
-    expr *child = (expr *) PyCapsule_GetPointer(child_capsule, EXPR_CAPSULE_NAME);
-    if (!child)
-    {
-        PyErr_SetString(PyExc_ValueError, "invalid child capsule");
-        return NULL;
-    }
-
-    PyArrayObject *a_array =
-        (PyArrayObject *) PyArray_FROM_OTF(a_obj, NPY_DOUBLE, NPY_ARRAY_IN_ARRAY);
-
-    if (!a_array)
-    {
-        return NULL;
-    }
-
-    /* Verify size matches child size */
-    int a_size = (int) PyArray_SIZE(a_array);
-    if (a_size != child->size)
-    {
-        Py_DECREF(a_array);
-        PyErr_SetString(PyExc_ValueError, "vector 'a' size must match child size");
-        return NULL;
-    }
-
-    double *a_data = (double *) PyArray_DATA(a_array);
-
-    expr *node = new_const_vector_mult(a_data, child);
-
-    Py_DECREF(a_array);
-
-    if (!node)
-    {
-        PyErr_SetString(PyExc_RuntimeError,
-                        "failed to create const_vector_mult node");
-        return NULL;
-    }
-    expr_retain(node); /* Capsule owns a reference */
-    return PyCapsule_New(node, EXPR_CAPSULE_NAME, expr_capsule_destructor);
-}
-
-#endif /* ATOM_CONST_VECTOR_MULT_H */
diff --git a/python/atoms/constant.h b/python/atoms/constant.h
deleted file mode 100644
index d5fcba3..0000000
--- a/python/atoms/constant.h
+++ /dev/null
@@ -1,35 +0,0 @@
-#ifndef ATOM_CONSTANT_H
-#define ATOM_CONSTANT_H
-
-#include "common.h"
-
-static PyObject *py_make_constant(PyObject *self, PyObject *args)
-{
-    int d1, d2, n_vars;
-    PyObject *values_obj;
-    if (!PyArg_ParseTuple(args, "iiiO", &d1, &d2, &n_vars, &values_obj))
-    {
-        return NULL;
-    }
-
-    PyArrayObject *values_array = (PyArrayObject *) PyArray_FROM_OTF(
-        values_obj, NPY_DOUBLE, NPY_ARRAY_IN_ARRAY);
-    if (!values_array)
-    {
-        return NULL;
-    }
-
-    expr *node =
-        new_constant(d1, d2, n_vars, (const double *) PyArray_DATA(values_array));
-    Py_DECREF(values_array);
-
-    if (!node)
-    {
-        PyErr_SetString(PyExc_RuntimeError, "failed to create constant node");
-        return NULL;
-    }
-    expr_retain(node); /* Capsule owns a reference */
-    return PyCapsule_New(node, EXPR_CAPSULE_NAME, expr_capsule_destructor);
-}
-
-#endif /* ATOM_CONSTANT_H */
diff --git a/python/atoms/cos.h b/python/atoms/cos.h
deleted file mode 100644
index 3f798be..0000000
--- a/python/atoms/cos.h
+++ /dev/null
@@ -1,31 +0,0 @@
-#ifndef ATOM_COS_H
-#define ATOM_COS_H
-
-#include "common.h"
-#include "elementwise_univariate.h"
-
-static PyObject *py_make_cos(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_cos(child);
-    if (!node)
-    {
-        PyErr_SetString(PyExc_RuntimeError, "failed to create cos node");
-        return NULL;
-    }
-    expr_retain(node); /* Capsule owns a reference */
-    return PyCapsule_New(node, EXPR_CAPSULE_NAME, expr_capsule_destructor);
-}
-
-#endif /* ATOM_COS_H */
diff --git a/python/atoms/diag_vec.h b/python/atoms/diag_vec.h
deleted file mode 100644
index b1056d6..0000000
--- a/python/atoms/diag_vec.h
+++ /dev/null
@@ -1,32 +0,0 @@
-#ifndef ATOM_DIAG_VEC_H
-#define ATOM_DIAG_VEC_H
-
-#include "common.h"
-
-static PyObject *py_make_diag_vec(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)
-    {
-        return NULL;
-    }
-
-    expr *node = new_diag_vec(child);
-    if (!node)
-    {
-        PyErr_SetString(PyExc_RuntimeError, "failed to create diag_vec node");
-        return NULL;
-    }
-
-    expr_retain(node);
-    return PyCapsule_New(node, EXPR_CAPSULE_NAME, expr_capsule_destructor);
-}
-
-#endif /* ATOM_DIAG_VEC_H */
diff --git a/python/atoms/entr.h b/python/atoms/entr.h
deleted file mode 100644
index 9a37668..0000000
--- a/python/atoms/entr.h
+++ /dev/null
@@ -1,31 +0,0 @@
-#ifndef ATOM_ENTR_H
-#define ATOM_ENTR_H
-
-#include "common.h"
-#include "elementwise_univariate.h"
-
-static PyObject *py_make_entr(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_entr(child);
-    if (!node)
-    {
-        PyErr_SetString(PyExc_RuntimeError, "failed to create entr node");
-        return NULL;
-    }
-    expr_retain(node); /* Capsule owns a reference */
-    return PyCapsule_New(node, EXPR_CAPSULE_NAME, expr_capsule_destructor);
-}
-
-#endif /* ATOM_ENTR_H */
diff --git a/python/atoms/exp.h b/python/atoms/exp.h
deleted file mode 100644
index 9ac7fc5..0000000
--- a/python/atoms/exp.h
+++ /dev/null
@@ -1,30 +0,0 @@
-#ifndef ATOM_EXP_H
-#define ATOM_EXP_H
-
-#include "common.h"
-
-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;
-    }
-    expr_retain(node); /* Capsule owns a reference */
-    return PyCapsule_New(node, EXPR_CAPSULE_NAME, expr_capsule_destructor);
-}
-
-#endif /* ATOM_EXP_H */
diff --git a/python/atoms/getters.h b/python/atoms/getters.h
deleted file mode 100644
index 538680a..0000000
--- a/python/atoms/getters.h
+++ /dev/null
@@ -1,43 +0,0 @@
-#ifndef ATOM_GETTERS_H
-#define ATOM_GETTERS_H
-
-#include "common.h"
-
-static PyObject *py_get_expr_dimensions(PyObject *self, PyObject *args)
-{
-    PyObject *expr_capsule;
-
-    if (!PyArg_ParseTuple(args, "O", &expr_capsule))
-    {
-        return NULL;
-    }
-
-    expr *node = (expr *) PyCapsule_GetPointer(expr_capsule, EXPR_CAPSULE_NAME);
-    if (!node)
-    {
-        return NULL;
-    }
-
-    // Return tuple (d1, d2)
-    return Py_BuildValue("(ii)", node->d1, node->d2);
-}
-
-static PyObject *py_get_expr_size(PyObject *self, PyObject *args)
-{
-    PyObject *expr_capsule;
-
-    if (!PyArg_ParseTuple(args, "O", &expr_capsule))
-    {
-        return NULL;
-    }
-
-    expr *node = (expr *) PyCapsule_GetPointer(expr_capsule, EXPR_CAPSULE_NAME);
-    if (!node)
-    {
-        return NULL;
-    }
-
-    return Py_BuildValue("i", node->size);
-}
-
-#endif /* ATOM_GETTERS_H */
diff --git a/python/atoms/hstack.h b/python/atoms/hstack.h
deleted file mode 100644
index e1ac188..0000000
--- a/python/atoms/hstack.h
+++ /dev/null
@@ -1,50 +0,0 @@
-#ifndef ATOM_HSTACK_H
-#define ATOM_HSTACK_H
-
-#include "common.h"
-
-static PyObject *py_make_hstack(PyObject *self, PyObject *args)
-{
-    PyObject *list_obj;
-    if (!PyArg_ParseTuple(args, "O", &list_obj))
-    {
-        return NULL;
-    }
-    if (!PyList_Check(list_obj))
-    {
-        PyErr_SetString(PyExc_TypeError,
-                        "First argument must be a list of expr capsules");
-        return NULL;
-    }
-    Py_ssize_t n_args = PyList_Size(list_obj);
-    if (n_args == 0)
-    {
-        PyErr_SetString(PyExc_ValueError, "List of expr capsules cannot be empty");
-        return NULL;
-    }
-    expr **expr_args = (expr **) calloc(n_args, sizeof(expr *));
-    for (Py_ssize_t i = 0; i < n_args; ++i)
-    {
-        PyObject *item = PyList_GetItem(list_obj, i);
-        expr *e = (expr *) PyCapsule_GetPointer(item, EXPR_CAPSULE_NAME);
-        if (!e)
-        {
-            free(expr_args);
-            PyErr_SetString(PyExc_ValueError, "Invalid expr capsule in list");
-            return NULL;
-        }
-        expr_args[i] = e;
-    }
-    int n_vars = expr_args[0]->n_vars;
-    expr *node = new_hstack(expr_args, (int) n_args, n_vars);
-    free(expr_args);
-    if (!node)
-    {
-        PyErr_SetString(PyExc_RuntimeError, "failed to create hstack node");
-        return NULL;
-    }
-    expr_retain(node);
-    return PyCapsule_New(node, EXPR_CAPSULE_NAME, expr_capsule_destructor);
-}
-
-#endif // ATOM_HSTACK_H
diff --git a/python/atoms/index.h b/python/atoms/index.h
deleted file mode 100644
index 64125ad..0000000
--- a/python/atoms/index.h
+++ /dev/null
@@ -1,54 +0,0 @@
-// SPDX-License-Identifier: Apache-2.0
-
-#ifndef ATOM_INDEX_H
-#define ATOM_INDEX_H
-
-#include "affine.h"
-#include "common.h"
-
-/* Index/slicing: y = child[indices] where indices is a list of flattened positions
- */
-static PyObject *py_make_index(PyObject *self, PyObject *args)
-{
-    PyObject *child_capsule;
-    int d1, d2;
-    PyObject *indices_obj;
-
-    if (!PyArg_ParseTuple(args, "OiiO", &child_capsule, &d1, &d2, &indices_obj))
-    {
-        return NULL;
-    }
-
-    expr *child = (expr *) PyCapsule_GetPointer(child_capsule, EXPR_CAPSULE_NAME);
-    if (!child)
-    {
-        PyErr_SetString(PyExc_ValueError, "invalid child capsule");
-        return NULL;
-    }
-
-    /* Convert indices array to int32 */
-    PyArrayObject *indices_array = (PyArrayObject *) PyArray_FROM_OTF(
-        indices_obj, NPY_INT32, NPY_ARRAY_IN_ARRAY);
-
-    if (!indices_array)
-    {
-        return NULL;
-    }
-
-    int n_idxs = (int) PyArray_SIZE(indices_array);
-    int *indices_data = (int *) PyArray_DATA(indices_array);
-
-    expr *node = new_index(child, d1, d2, indices_data, n_idxs);
-
-    Py_DECREF(indices_array);
-
-    if (!node)
-    {
-        PyErr_SetString(PyExc_RuntimeError, "failed to create index node");
-        return NULL;
-    }
-    expr_retain(node); /* Capsule owns a reference */
-    return PyCapsule_New(node, EXPR_CAPSULE_NAME, expr_capsule_destructor);
-}
-
-#endif /* ATOM_INDEX_H */
diff --git a/python/atoms/left_matmul.h b/python/atoms/left_matmul.h
deleted file mode 100644
index 27fe3a4..0000000
--- a/python/atoms/left_matmul.h
+++ /dev/null
@@ -1,63 +0,0 @@
-#ifndef ATOM_LEFT_MATMUL_H
-#define ATOM_LEFT_MATMUL_H
-
-#include "bivariate.h"
-#include "common.h"
-
-/* Left matrix multiplication: A @ f(x) where A is a constant matrix */
-static PyObject *py_make_left_matmul(PyObject *self, PyObject *args)
-{
-    PyObject *child_capsule;
-    PyObject *data_obj, *indices_obj, *indptr_obj;
-    int m, n;
-    if (!PyArg_ParseTuple(args, "OOOOii", &child_capsule, &data_obj, &indices_obj,
-                          &indptr_obj, &m, &n))
-    {
-        return NULL;
-    }
-
-    expr *child = (expr *) PyCapsule_GetPointer(child_capsule, EXPR_CAPSULE_NAME);
-    if (!child)
-    {
-        PyErr_SetString(PyExc_ValueError, "invalid child capsule");
-        return NULL;
-    }
-
-    PyArrayObject *data_array =
-        (PyArrayObject *) PyArray_FROM_OTF(data_obj, NPY_DOUBLE, NPY_ARRAY_IN_ARRAY);
-    PyArrayObject *indices_array = (PyArrayObject *) PyArray_FROM_OTF(
-        indices_obj, NPY_INT32, NPY_ARRAY_IN_ARRAY);
-    PyArrayObject *indptr_array = (PyArrayObject *) PyArray_FROM_OTF(
-        indptr_obj, NPY_INT32, NPY_ARRAY_IN_ARRAY);
-
-    if (!data_array || !indices_array || !indptr_array)
-    {
-        Py_XDECREF(data_array);
-        Py_XDECREF(indices_array);
-        Py_XDECREF(indptr_array);
-        return NULL;
-    }
-
-    int nnz = (int) PyArray_SIZE(data_array);
-    CSR_Matrix *A = new_csr_matrix(m, n, nnz);
-    memcpy(A->x, PyArray_DATA(data_array), nnz * sizeof(double));
-    memcpy(A->i, PyArray_DATA(indices_array), nnz * sizeof(int));
-    memcpy(A->p, PyArray_DATA(indptr_array), (m + 1) * sizeof(int));
-
-    Py_DECREF(data_array);
-    Py_DECREF(indices_array);
-    Py_DECREF(indptr_array);
-
-    expr *node = new_left_matmul(child, A);
-    free_csr_matrix(A);
-
-    if (!node)
-    {
-        PyErr_SetString(PyExc_RuntimeError, "failed to create left_matmul node");
-        return NULL;
-    }
-    expr_retain(node); /* Capsule owns a reference */
-    return PyCapsule_New(node, EXPR_CAPSULE_NAME, expr_capsule_destructor);
-}
-
-#endif /* ATOM_LEFT_MATMUL_H */
diff --git a/python/atoms/linear.h b/python/atoms/linear.h
deleted file mode 100644
index 3465e05..0000000
--- a/python/atoms/linear.h
+++ /dev/null
@@ -1,83 +0,0 @@
-#ifndef ATOM_LINEAR_H
-#define ATOM_LINEAR_H
-
-#include "common.h"
-
-static PyObject *py_make_linear(PyObject *self, PyObject *args)
-{
-    PyObject *child_capsule;
-    PyObject *data_obj, *indices_obj, *indptr_obj;
-    PyObject *b_obj = Py_None; /* Optional offset vector */
-    int m, n;
-
-    /* Accept optional b array: (child, data, indices, indptr, m, n[, b]) */
-    if (!PyArg_ParseTuple(args, "OOOOii|O", &child_capsule, &data_obj, &indices_obj,
-                          &indptr_obj, &m, &n, &b_obj))
-    {
-        return NULL;
-    }
-
-    expr *child = (expr *) PyCapsule_GetPointer(child_capsule, EXPR_CAPSULE_NAME);
-    if (!child)
-    {
-        PyErr_SetString(PyExc_ValueError, "invalid child capsule");
-        return NULL;
-    }
-
-    PyArrayObject *data_array =
-        (PyArrayObject *) PyArray_FROM_OTF(data_obj, NPY_DOUBLE, NPY_ARRAY_IN_ARRAY);
-    PyArrayObject *indices_array = (PyArrayObject *) PyArray_FROM_OTF(
-        indices_obj, NPY_INT32, NPY_ARRAY_IN_ARRAY);
-    PyArrayObject *indptr_array = (PyArrayObject *) PyArray_FROM_OTF(
-        indptr_obj, NPY_INT32, NPY_ARRAY_IN_ARRAY);
-
-    if (!data_array || !indices_array || !indptr_array)
-    {
-        Py_XDECREF(data_array);
-        Py_XDECREF(indices_array);
-        Py_XDECREF(indptr_array);
-        return NULL;
-    }
-
-    int nnz = (int) PyArray_SIZE(data_array);
-    CSR_Matrix *A = new_csr_matrix(m, n, nnz);
-    memcpy(A->x, PyArray_DATA(data_array), nnz * sizeof(double));
-    memcpy(A->i, PyArray_DATA(indices_array), nnz * sizeof(int));
-    memcpy(A->p, PyArray_DATA(indptr_array), (m + 1) * sizeof(int));
-
-    Py_DECREF(data_array);
-    Py_DECREF(indices_array);
-    Py_DECREF(indptr_array);
-
-    /* Parse optional b offset vector */
-    double *b_data = NULL;
-    PyArrayObject *b_array = NULL;
-
-    if (b_obj != Py_None)
-    {
-        b_array = (PyArrayObject *) PyArray_FROM_OTF(b_obj, NPY_DOUBLE,
-                                                     NPY_ARRAY_IN_ARRAY);
-        if (!b_array)
-        {
-            free_csr_matrix(A);
-            return NULL;
-        }
-        b_data = (double *) PyArray_DATA(b_array);
-    }
-
-    expr *node = new_linear(child, A, b_data);
-
-    /* Clean up */
-    free_csr_matrix(A);
-    Py_XDECREF(b_array);
-
-    if (!node)
-    {
-        PyErr_SetString(PyExc_RuntimeError, "failed to create linear node");
-        return NULL;
-    }
-    expr_retain(node); /* Capsule owns a reference */
-    return PyCapsule_New(node, EXPR_CAPSULE_NAME, expr_capsule_destructor);
-}
-
-#endif /* ATOM_LINEAR_H */
diff --git a/python/atoms/log.h b/python/atoms/log.h
deleted file mode 100644
index 738b4e5..0000000
--- a/python/atoms/log.h
+++ /dev/null
@@ -1,30 +0,0 @@
-#ifndef ATOM_LOG_H
-#define ATOM_LOG_H
-
-#include "common.h"
-
-static PyObject *py_make_log(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_log(child);
-    if (!node)
-    {
-        PyErr_SetString(PyExc_RuntimeError, "failed to create log node");
-        return NULL;
-    }
-    expr_retain(node); /* Capsule owns a reference */
-    return PyCapsule_New(node, EXPR_CAPSULE_NAME, expr_capsule_destructor);
-}
-
-#endif /* ATOM_LOG_H */
diff --git a/python/atoms/logistic.h b/python/atoms/logistic.h
deleted file mode 100644
index e903661..0000000
--- a/python/atoms/logistic.h
+++ /dev/null
@@ -1,31 +0,0 @@
-#ifndef ATOM_LOGISTIC_H
-#define ATOM_LOGISTIC_H
-
-#include "common.h"
-#include "elementwise_univariate.h"
-
-static PyObject *py_make_logistic(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_logistic(child);
-    if (!node)
-    {
-        PyErr_SetString(PyExc_RuntimeError, "failed to create logistic node");
-        return NULL;
-    }
-    expr_retain(node); /* Capsule owns a reference */
-    return PyCapsule_New(node, EXPR_CAPSULE_NAME, expr_capsule_destructor);
-}
-
-#endif /* ATOM_LOGISTIC_H */
diff --git a/python/atoms/matmul.h b/python/atoms/matmul.h
deleted file mode 100644
index 9bfe859..0000000
--- a/python/atoms/matmul.h
+++ /dev/null
@@ -1,33 +0,0 @@
-#ifndef ATOM_MATMUL_H
-#define ATOM_MATMUL_H
-
-#include "bivariate.h"
-#include "common.h"
-
-/* Matrix multiplication: Z = X @ Y */
-static PyObject *py_make_matmul(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_matmul(left, right);
-    if (!node)
-    {
-        PyErr_SetString(PyExc_RuntimeError, "failed to create matmul node");
-        return NULL;
-    }
-    expr_retain(node); /* Capsule owns a reference */
-    return PyCapsule_New(node, EXPR_CAPSULE_NAME, expr_capsule_destructor);
-}
-
-#endif /* ATOM_MATMUL_H */
diff --git a/python/atoms/multiply.h b/python/atoms/multiply.h
deleted file mode 100644
index 5403b38..0000000
--- a/python/atoms/multiply.h
+++ /dev/null
@@ -1,32 +0,0 @@
-#ifndef ATOM_MULTIPLY_H
-#define ATOM_MULTIPLY_H
-
-#include "bivariate.h"
-#include "common.h"
-
-static PyObject *py_make_multiply(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_elementwise_mult(left, right);
-    if (!node)
-    {
-        PyErr_SetString(PyExc_RuntimeError, "failed to create multiply node");
-        return NULL;
-    }
-    expr_retain(node); /* Capsule owns a reference */
-    return PyCapsule_New(node, EXPR_CAPSULE_NAME, expr_capsule_destructor);
-}
-
-#endif /* ATOM_MULTIPLY_H */
diff --git a/python/atoms/neg.h b/python/atoms/neg.h
deleted file mode 100644
index 3927b0c..0000000
--- a/python/atoms/neg.h
+++ /dev/null
@@ -1,30 +0,0 @@
-#ifndef ATOM_NEG_H
-#define ATOM_NEG_H
-
-#include "common.h"
-
-static PyObject *py_make_neg(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_neg(child);
-    if (!node)
-    {
-        PyErr_SetString(PyExc_RuntimeError, "failed to create neg node");
-        return NULL;
-    }
-    expr_retain(node); /* Capsule owns a reference */
-    return PyCapsule_New(node, EXPR_CAPSULE_NAME, expr_capsule_destructor);
-}
-
-#endif /* ATOM_NEG_H */
diff --git a/python/atoms/power.h b/python/atoms/power.h
deleted file mode 100644
index da9b7bf..0000000
--- a/python/atoms/power.h
+++ /dev/null
@@ -1,32 +0,0 @@
-#ifndef ATOM_POWER_H
-#define ATOM_POWER_H
-
-#include "common.h"
-#include "elementwise_univariate.h"
-
-static PyObject *py_make_power(PyObject *self, PyObject *args)
-{
-    PyObject *child_capsule;
-    double p;
-    if (!PyArg_ParseTuple(args, "Od", &child_capsule, &p))
-    {
-        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_power(child, p);
-    if (!node)
-    {
-        PyErr_SetString(PyExc_RuntimeError, "failed to create power node");
-        return NULL;
-    }
-    expr_retain(node); /* Capsule owns a reference */
-    return PyCapsule_New(node, EXPR_CAPSULE_NAME, expr_capsule_destructor);
-}
-
-#endif /* ATOM_POWER_H */
diff --git a/python/atoms/prod.h b/python/atoms/prod.h
deleted file mode 100644
index 9e15818..0000000
--- a/python/atoms/prod.h
+++ /dev/null
@@ -1,33 +0,0 @@
-// SPDX-License-Identifier: Apache-2.0
-#ifndef ATOM_PROD_H
-#define ATOM_PROD_H
-
-#include "common.h"
-#include "other.h"
-
-static PyObject *py_make_prod(PyObject *self, PyObject *args)
-{
-    (void) self;
-    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_prod(child);
-    if (!node)
-    {
-        PyErr_SetString(PyExc_RuntimeError, "failed to create prod node");
-        return NULL;
-    }
-    expr_retain(node); /* Capsule owns a reference */
-    return PyCapsule_New(node, EXPR_CAPSULE_NAME, expr_capsule_destructor);
-}
-
-#endif /* ATOM_PROD_H */
diff --git a/python/atoms/prod_axis_one.h b/python/atoms/prod_axis_one.h
deleted file mode 100644
index d9c3b8a..0000000
--- a/python/atoms/prod_axis_one.h
+++ /dev/null
@@ -1,32 +0,0 @@
-#ifndef ATOM_PROD_AXIS_ONE_H
-#define ATOM_PROD_AXIS_ONE_H
-
-#include "common.h"
-#include "other.h"
-
-static PyObject *py_make_prod_axis_one(PyObject *self, PyObject *args)
-{
-    (void) self;
-    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_prod_axis_one(child);
-    if (!node)
-    {
-        PyErr_SetString(PyExc_RuntimeError, "failed to create prod_axis_one node");
-        return NULL;
-    }
-    expr_retain(node); /* Capsule owns a reference */
-    return PyCapsule_New(node, EXPR_CAPSULE_NAME, expr_capsule_destructor);
-}
-
-#endif /* ATOM_PROD_AXIS_ONE_H */
diff --git a/python/atoms/prod_axis_zero.h b/python/atoms/prod_axis_zero.h
deleted file mode 100644
index 02b30f0..0000000
--- a/python/atoms/prod_axis_zero.h
+++ /dev/null
@@ -1,33 +0,0 @@
-// SPDX-License-Identifier: Apache-2.0
-#ifndef ATOM_PROD_AXIS_ZERO_H
-#define ATOM_PROD_AXIS_ZERO_H
-
-#include "common.h"
-#include "other.h"
-
-static PyObject *py_make_prod_axis_zero(PyObject *self, PyObject *args)
-{
-    (void) self;
-    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_prod_axis_zero(child);
-    if (!node)
-    {
-        PyErr_SetString(PyExc_RuntimeError, "failed to create prod_axis_zero node");
-        return NULL;
-    }
-    expr_retain(node); /* Capsule owns a reference */
-    return PyCapsule_New(node, EXPR_CAPSULE_NAME, expr_capsule_destructor);
-}
-
-#endif /* ATOM_PROD_AXIS_ZERO_H */
diff --git a/python/atoms/promote.h b/python/atoms/promote.h
deleted file mode 100644
index 97e405c..0000000
--- a/python/atoms/promote.h
+++ /dev/null
@@ -1,31 +0,0 @@
-#ifndef ATOM_PROMOTE_H
-#define ATOM_PROMOTE_H
-
-#include "common.h"
-
-static PyObject *py_make_promote(PyObject *self, PyObject *args)
-{
-    PyObject *child_capsule;
-    int d1, d2;
-    if (!PyArg_ParseTuple(args, "Oii", &child_capsule, &d1, &d2))
-    {
-        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_promote(child, d1, d2);
-    if (!node)
-    {
-        PyErr_SetString(PyExc_RuntimeError, "failed to create promote node");
-        return NULL;
-    }
-    expr_retain(node); /* Capsule owns a reference */
-    return PyCapsule_New(node, EXPR_CAPSULE_NAME, expr_capsule_destructor);
-}
-
-#endif /* ATOM_PROMOTE_H */
diff --git a/python/atoms/quad_form.h b/python/atoms/quad_form.h
deleted file mode 100644
index 036640a..0000000
--- a/python/atoms/quad_form.h
+++ /dev/null
@@ -1,64 +0,0 @@
-#ifndef ATOM_QUAD_FORM_H
-#define ATOM_QUAD_FORM_H
-
-#include "common.h"
-#include "other.h"
-
-/* Quadratic form: y = x' * Q * x where Q is a constant matrix */
-static PyObject *py_make_quad_form(PyObject *self, PyObject *args)
-{
-    PyObject *child_capsule;
-    PyObject *data_obj, *indices_obj, *indptr_obj;
-    int m, n;
-
-    if (!PyArg_ParseTuple(args, "OOOOii", &child_capsule, &data_obj, &indices_obj,
-                          &indptr_obj, &m, &n))
-    {
-        return NULL;
-    }
-
-    expr *child = (expr *) PyCapsule_GetPointer(child_capsule, EXPR_CAPSULE_NAME);
-    if (!child)
-    {
-        PyErr_SetString(PyExc_ValueError, "invalid child capsule");
-        return NULL;
-    }
-
-    PyArrayObject *data_array =
-        (PyArrayObject *) PyArray_FROM_OTF(data_obj, NPY_DOUBLE, NPY_ARRAY_IN_ARRAY);
-    PyArrayObject *indices_array = (PyArrayObject *) PyArray_FROM_OTF(
-        indices_obj, NPY_INT32, NPY_ARRAY_IN_ARRAY);
-    PyArrayObject *indptr_array = (PyArrayObject *) PyArray_FROM_OTF(
-        indptr_obj, NPY_INT32, NPY_ARRAY_IN_ARRAY);
-
-    if (!data_array || !indices_array || !indptr_array)
-    {
-        Py_XDECREF(data_array);
-        Py_XDECREF(indices_array);
-        Py_XDECREF(indptr_array);
-        return NULL;
-    }
-
-    int nnz = (int) PyArray_SIZE(data_array);
-    CSR_Matrix *Q = new_csr_matrix(m, n, nnz);
-    memcpy(Q->x, PyArray_DATA(data_array), nnz * sizeof(double));
-    memcpy(Q->i, PyArray_DATA(indices_array), nnz * sizeof(int));
-    memcpy(Q->p, PyArray_DATA(indptr_array), (m + 1) * sizeof(int));
-
-    Py_DECREF(data_array);
-    Py_DECREF(indices_array);
-    Py_DECREF(indptr_array);
-
-    expr *node = new_quad_form(child, Q);
-    free_csr_matrix(Q);
-
-    if (!node)
-    {
-        PyErr_SetString(PyExc_RuntimeError, "failed to create quad_form node");
-        return NULL;
-    }
-    expr_retain(node); /* Capsule owns a reference */
-    return PyCapsule_New(node, EXPR_CAPSULE_NAME, expr_capsule_destructor);
-}
-
-#endif /* ATOM_QUAD_FORM_H */
diff --git a/python/atoms/quad_over_lin.h b/python/atoms/quad_over_lin.h
deleted file mode 100644
index d063c8a..0000000
--- a/python/atoms/quad_over_lin.h
+++ /dev/null
@@ -1,36 +0,0 @@
-// SPDX-License-Identifier: Apache-2.0
-
-#ifndef ATOM_QUAD_OVER_LIN_H
-#define ATOM_QUAD_OVER_LIN_H
-
-#include "bivariate.h"
-#include "common.h"
-
-/* quad_over_lin: y = sum(x^2) / z where x is left, z is right (scalar) */
-static PyObject *py_make_quad_over_lin(PyObject *self, PyObject *args)
-{
-    (void) self;
-    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_quad_over_lin(left, right);
-    if (!node)
-    {
-        PyErr_SetString(PyExc_RuntimeError, "failed to create quad_over_lin node");
-        return NULL;
-    }
-    expr_retain(node); /* Capsule owns a reference */
-    return PyCapsule_New(node, EXPR_CAPSULE_NAME, expr_capsule_destructor);
-}
-
-#endif /* ATOM_QUAD_OVER_LIN_H */
diff --git a/python/atoms/rel_entr.h b/python/atoms/rel_entr.h
deleted file mode 100644
index 1b74318..0000000
--- a/python/atoms/rel_entr.h
+++ /dev/null
@@ -1,36 +0,0 @@
-// SPDX-License-Identifier: Apache-2.0
-
-#ifndef ATOM_REL_ENTR_H
-#define ATOM_REL_ENTR_H
-
-#include "bivariate.h"
-#include "common.h"
-
-/* rel_entr: rel_entr(x, y) = x * log(x/y) elementwise */
-static PyObject *py_make_rel_entr(PyObject *self, PyObject *args)
-{
-    (void) self;
-    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_rel_entr_vector_args(left, right);
-    if (!node)
-    {
-        PyErr_SetString(PyExc_RuntimeError, "failed to create rel_entr node");
-        return NULL;
-    }
-    expr_retain(node); /* Capsule owns a reference */
-    return PyCapsule_New(node, EXPR_CAPSULE_NAME, expr_capsule_destructor);
-}
-
-#endif /* ATOM_REL_ENTR_H */
diff --git a/python/atoms/rel_entr_scalar_vector.h b/python/atoms/rel_entr_scalar_vector.h
deleted file mode 100644
index 85bcfe7..0000000
--- a/python/atoms/rel_entr_scalar_vector.h
+++ /dev/null
@@ -1,37 +0,0 @@
-// SPDX-License-Identifier: Apache-2.0
-
-#ifndef ATOM_REL_ENTR_SCALAR_VECTOR_H
-#define ATOM_REL_ENTR_SCALAR_VECTOR_H
-
-#include "bivariate.h"
-#include "common.h"
-
-/* rel_entr_scalar_vector: rel_entr(x, y) where x is scalar, y is vector */
-static PyObject *py_make_rel_entr_scalar_vector(PyObject *self, PyObject *args)
-{
-    (void) self;
-    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_rel_entr_first_arg_scalar(left, right);
-    if (!node)
-    {
-        PyErr_SetString(PyExc_RuntimeError,
-                        "failed to create rel_entr_scalar_vector node");
-        return NULL;
-    }
-    expr_retain(node); /* Capsule owns a reference */
-    return PyCapsule_New(node, EXPR_CAPSULE_NAME, expr_capsule_destructor);
-}
-
-#endif /* ATOM_REL_ENTR_SCALAR_VECTOR_H */
diff --git a/python/atoms/rel_entr_vector_scalar.h b/python/atoms/rel_entr_vector_scalar.h
deleted file mode 100644
index d671646..0000000
--- a/python/atoms/rel_entr_vector_scalar.h
+++ /dev/null
@@ -1,37 +0,0 @@
-// SPDX-License-Identifier: Apache-2.0
-
-#ifndef ATOM_REL_ENTR_VECTOR_SCALAR_H
-#define ATOM_REL_ENTR_VECTOR_SCALAR_H
-
-#include "bivariate.h"
-#include "common.h"
-
-/* rel_entr_vector_scalar: rel_entr(x, y) where x is vector, y is scalar */
-static PyObject *py_make_rel_entr_vector_scalar(PyObject *self, PyObject *args)
-{
-    (void) self;
-    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_rel_entr_second_arg_scalar(left, right);
-    if (!node)
-    {
-        PyErr_SetString(PyExc_RuntimeError,
-                        "failed to create rel_entr_vector_scalar node");
-        return NULL;
-    }
-    expr_retain(node); /* Capsule owns a reference */
-    return PyCapsule_New(node, EXPR_CAPSULE_NAME, expr_capsule_destructor);
-}
-
-#endif /* ATOM_REL_ENTR_VECTOR_SCALAR_H */
diff --git a/python/atoms/reshape.h b/python/atoms/reshape.h
deleted file mode 100644
index 8547a79..0000000
--- a/python/atoms/reshape.h
+++ /dev/null
@@ -1,33 +0,0 @@
-#ifndef ATOM_RESHAPE_H
-#define ATOM_RESHAPE_H
-
-#include "common.h"
-
-static PyObject *py_make_reshape(PyObject *self, PyObject *args)
-{
-    PyObject *child_capsule;
-    int d1, d2;
-
-    if (!PyArg_ParseTuple(args, "Oii", &child_capsule, &d1, &d2))
-    {
-        return NULL;
-    }
-
-    expr *child = (expr *) PyCapsule_GetPointer(child_capsule, EXPR_CAPSULE_NAME);
-    if (!child)
-    {
-        return NULL;
-    }
-
-    expr *node = new_reshape(child, d1, d2);
-    if (!node)
-    {
-        PyErr_SetString(PyExc_RuntimeError, "failed to create reshape node");
-        return NULL;
-    }
-
-    expr_retain(node);
-    return PyCapsule_New(node, EXPR_CAPSULE_NAME, expr_capsule_destructor);
-}
-
-#endif /* ATOM_RESHAPE_H */
diff --git a/python/atoms/right_matmul.h b/python/atoms/right_matmul.h
deleted file mode 100644
index c1c3481..0000000
--- a/python/atoms/right_matmul.h
+++ /dev/null
@@ -1,63 +0,0 @@
-#ifndef ATOM_RIGHT_MATMUL_H
-#define ATOM_RIGHT_MATMUL_H
-
-#include "bivariate.h"
-#include "common.h"
-
-/* Right matrix multiplication: f(x) @ A where A is a constant matrix */
-static PyObject *py_make_right_matmul(PyObject *self, PyObject *args)
-{
-    PyObject *child_capsule;
-    PyObject *data_obj, *indices_obj, *indptr_obj;
-    int m, n;
-    if (!PyArg_ParseTuple(args, "OOOOii", &child_capsule, &data_obj, &indices_obj,
-                          &indptr_obj, &m, &n))
-    {
-        return NULL;
-    }
-
-    expr *child = (expr *) PyCapsule_GetPointer(child_capsule, EXPR_CAPSULE_NAME);
-    if (!child)
-    {
-        PyErr_SetString(PyExc_ValueError, "invalid child capsule");
-        return NULL;
-    }
-
-    PyArrayObject *data_array =
-        (PyArrayObject *) PyArray_FROM_OTF(data_obj, NPY_DOUBLE, NPY_ARRAY_IN_ARRAY);
-    PyArrayObject *indices_array = (PyArrayObject *) PyArray_FROM_OTF(
-        indices_obj, NPY_INT32, NPY_ARRAY_IN_ARRAY);
-    PyArrayObject *indptr_array = (PyArrayObject *) PyArray_FROM_OTF(
-        indptr_obj, NPY_INT32, NPY_ARRAY_IN_ARRAY);
-
-    if (!data_array || !indices_array || !indptr_array)
-    {
-        Py_XDECREF(data_array);
-        Py_XDECREF(indices_array);
-        Py_XDECREF(indptr_array);
-        return NULL;
-    }
-
-    int nnz = (int) PyArray_SIZE(data_array);
-    CSR_Matrix *A = new_csr_matrix(m, n, nnz);
-    memcpy(A->x, PyArray_DATA(data_array), nnz * sizeof(double));
-    memcpy(A->i, PyArray_DATA(indices_array), nnz * sizeof(int));
-    memcpy(A->p, PyArray_DATA(indptr_array), (m + 1) * sizeof(int));
-
-    Py_DECREF(data_array);
-    Py_DECREF(indices_array);
-    Py_DECREF(indptr_array);
-
-    expr *node = new_right_matmul(child, A);
-    free_csr_matrix(A);
-
-    if (!node)
-    {
-        PyErr_SetString(PyExc_RuntimeError, "failed to create right_matmul node");
-        return NULL;
-    }
-    expr_retain(node); /* Capsule owns a reference */
-    return PyCapsule_New(node, EXPR_CAPSULE_NAME, expr_capsule_destructor);
-}
-
-#endif /* ATOM_RIGHT_MATMUL_H */
diff --git a/python/atoms/sin.h b/python/atoms/sin.h
deleted file mode 100644
index 426d302..0000000
--- a/python/atoms/sin.h
+++ /dev/null
@@ -1,31 +0,0 @@
-#ifndef ATOM_SIN_H
-#define ATOM_SIN_H
-
-#include "common.h"
-#include "elementwise_univariate.h"
-
-static PyObject *py_make_sin(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_sin(child);
-    if (!node)
-    {
-        PyErr_SetString(PyExc_RuntimeError, "failed to create sin node");
-        return NULL;
-    }
-    expr_retain(node); /* Capsule owns a reference */
-    return PyCapsule_New(node, EXPR_CAPSULE_NAME, expr_capsule_destructor);
-}
-
-#endif /* ATOM_SIN_H */
diff --git a/python/atoms/sinh.h b/python/atoms/sinh.h
deleted file mode 100644
index be86cf5..0000000
--- a/python/atoms/sinh.h
+++ /dev/null
@@ -1,31 +0,0 @@
-#ifndef ATOM_SINH_H
-#define ATOM_SINH_H
-
-#include "common.h"
-#include "elementwise_univariate.h"
-
-static PyObject *py_make_sinh(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_sinh(child);
-    if (!node)
-    {
-        PyErr_SetString(PyExc_RuntimeError, "failed to create sinh node");
-        return NULL;
-    }
-    expr_retain(node); /* Capsule owns a reference */
-    return PyCapsule_New(node, EXPR_CAPSULE_NAME, expr_capsule_destructor);
-}
-
-#endif /* ATOM_SINH_H */
diff --git a/python/atoms/sum.h b/python/atoms/sum.h
deleted file mode 100644
index 8d01604..0000000
--- a/python/atoms/sum.h
+++ /dev/null
@@ -1,31 +0,0 @@
-#ifndef ATOM_SUM_H
-#define ATOM_SUM_H
-
-#include "common.h"
-
-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;
-    }
-    expr_retain(node); /* Capsule owns a reference */
-    return PyCapsule_New(node, EXPR_CAPSULE_NAME, expr_capsule_destructor);
-}
-
-#endif /* ATOM_SUM_H */
diff --git a/python/atoms/tan.h b/python/atoms/tan.h
deleted file mode 100644
index 3ec83b7..0000000
--- a/python/atoms/tan.h
+++ /dev/null
@@ -1,31 +0,0 @@
-#ifndef ATOM_TAN_H
-#define ATOM_TAN_H
-
-#include "common.h"
-#include "elementwise_univariate.h"
-
-static PyObject *py_make_tan(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_tan(child);
-    if (!node)
-    {
-        PyErr_SetString(PyExc_RuntimeError, "failed to create tan node");
-        return NULL;
-    }
-    expr_retain(node); /* Capsule owns a reference */
-    return PyCapsule_New(node, EXPR_CAPSULE_NAME, expr_capsule_destructor);
-}
-
-#endif /* ATOM_TAN_H */
diff --git a/python/atoms/tanh.h b/python/atoms/tanh.h
deleted file mode 100644
index e60d420..0000000
--- a/python/atoms/tanh.h
+++ /dev/null
@@ -1,31 +0,0 @@
-#ifndef ATOM_TANH_H
-#define ATOM_TANH_H
-
-#include "common.h"
-#include "elementwise_univariate.h"
-
-static PyObject *py_make_tanh(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_tanh(child);
-    if (!node)
-    {
-        PyErr_SetString(PyExc_RuntimeError, "failed to create tanh node");
-        return NULL;
-    }
-    expr_retain(node); /* Capsule owns a reference */
-    return PyCapsule_New(node, EXPR_CAPSULE_NAME, expr_capsule_destructor);
-}
-
-#endif /* ATOM_TANH_H */
diff --git a/python/atoms/trace.h b/python/atoms/trace.h
deleted file mode 100644
index e130aa2..0000000
--- a/python/atoms/trace.h
+++ /dev/null
@@ -1,29 +0,0 @@
-#ifndef ATOM_TRACE_H
-#define ATOM_TRACE_H
-
-#include "common.h"
-
-static PyObject *py_make_trace(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_trace(child);
-    if (!node)
-    {
-        PyErr_SetString(PyExc_RuntimeError, "failed to create trace node");
-        return NULL;
-    }
-    expr_retain(node); /* Capsule owns a reference */
-    return PyCapsule_New(node, EXPR_CAPSULE_NAME, expr_capsule_destructor);
-}
-
-#endif // ATOM_TRACE_H
diff --git a/python/atoms/transpose.h b/python/atoms/transpose.h
deleted file mode 100644
index 85c1dcc..0000000
--- a/python/atoms/transpose.h
+++ /dev/null
@@ -1,32 +0,0 @@
-#ifndef PYTHON_ATOMS_TRANSPOSE_H
-#define PYTHON_ATOMS_TRANSPOSE_H
-
-#include "common.h"
-
-// Python binding for the transpose atom
-static PyObject *py_make_transpose(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_transpose(child);
-    if (!node)
-    {
-        PyErr_SetString(PyExc_RuntimeError, "failed to create trace node");
-        return NULL;
-    }
-    expr_retain(node); /* Capsule owns a reference */
-    return PyCapsule_New(node, EXPR_CAPSULE_NAME, expr_capsule_destructor);
-}
-
-#endif // PYTHON_ATOMS_TRANSPOSE_H
diff --git a/python/atoms/variable.h b/python/atoms/variable.h
deleted file mode 100644
index 2dc1972..0000000
--- a/python/atoms/variable.h
+++ /dev/null
@@ -1,24 +0,0 @@
-#ifndef ATOM_VARIABLE_H
-#define ATOM_VARIABLE_H
-
-#include "common.h"
-
-static PyObject *py_make_variable(PyObject *self, PyObject *args)
-{
-    int d1, d2, var_id, n_vars;
-    if (!PyArg_ParseTuple(args, "iiii", &d1, &d2, &var_id, &n_vars))
-    {
-        return NULL;
-    }
-
-    expr *node = new_variable(d1, d2, var_id, n_vars);
-    if (!node)
-    {
-        PyErr_SetString(PyExc_RuntimeError, "failed to create variable node");
-        return NULL;
-    }
-    expr_retain(node); /* Capsule owns a reference */
-    return PyCapsule_New(node, EXPR_CAPSULE_NAME, expr_capsule_destructor);
-}
-
-#endif /* ATOM_VARIABLE_H */
diff --git a/python/atoms/xexp.h b/python/atoms/xexp.h
deleted file mode 100644
index 1dda312..0000000
--- a/python/atoms/xexp.h
+++ /dev/null
@@ -1,31 +0,0 @@
-#ifndef ATOM_XEXP_H
-#define ATOM_XEXP_H
-
-#include "common.h"
-#include "elementwise_univariate.h"
-
-static PyObject *py_make_xexp(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_xexp(child);
-    if (!node)
-    {
-        PyErr_SetString(PyExc_RuntimeError, "failed to create xexp node");
-        return NULL;
-    }
-    expr_retain(node); /* Capsule owns a reference */
-    return PyCapsule_New(node, EXPR_CAPSULE_NAME, expr_capsule_destructor);
-}
-
-#endif /* ATOM_XEXP_H */
diff --git a/python/bindings.c b/python/bindings.c
deleted file mode 100644
index d5b4713..0000000
--- a/python/bindings.c
+++ /dev/null
@@ -1,156 +0,0 @@
-#define PY_SSIZE_T_CLEAN
-#include <Python.h>
-#include <numpy/arrayobject.h>
-
-/* Include atom bindings */
-#include "atoms/add.h"
-#include "atoms/asinh.h"
-#include "atoms/atanh.h"
-#include "atoms/broadcast.h"
-#include "atoms/const_scalar_mult.h"
-#include "atoms/const_vector_mult.h"
-#include "atoms/constant.h"
-#include "atoms/cos.h"
-#include "atoms/diag_vec.h"
-#include "atoms/entr.h"
-#include "atoms/exp.h"
-#include "atoms/getters.h"
-#include "atoms/hstack.h"
-#include "atoms/index.h"
-#include "atoms/left_matmul.h"
-#include "atoms/linear.h"
-#include "atoms/log.h"
-#include "atoms/logistic.h"
-#include "atoms/matmul.h"
-#include "atoms/multiply.h"
-#include "atoms/neg.h"
-#include "atoms/power.h"
-#include "atoms/prod.h"
-#include "atoms/prod_axis_one.h"
-#include "atoms/prod_axis_zero.h"
-#include "atoms/promote.h"
-#include "atoms/quad_form.h"
-#include "atoms/quad_over_lin.h"
-#include "atoms/rel_entr.h"
-#include "atoms/rel_entr_scalar_vector.h"
-#include "atoms/rel_entr_vector_scalar.h"
-#include "atoms/reshape.h"
-#include "atoms/right_matmul.h"
-#include "atoms/sin.h"
-#include "atoms/sinh.h"
-#include "atoms/sum.h"
-#include "atoms/tan.h"
-#include "atoms/tanh.h"
-#include "atoms/trace.h"
-#include "atoms/transpose.h"
-#include "atoms/variable.h"
-#include "atoms/xexp.h"
-
-/* Include problem bindings */
-#include "problem/constraint_forward.h"
-#include "problem/gradient.h"
-#include "problem/hessian.h"
-#include "problem/init_derivatives.h"
-#include "problem/jacobian.h"
-#include "problem/make_problem.h"
-#include "problem/objective_forward.h"
-
-static int numpy_initialized = 0;
-
-static int ensure_numpy(void)
-{
-    if (numpy_initialized) return 0;
-    import_array1(-1);
-    numpy_initialized = 1;
-    return 0;
-}
-
-static PyMethodDef DNLPMethods[] = {
-    {"make_variable", py_make_variable, METH_VARARGS, "Create variable node"},
-    {"make_constant", py_make_constant, METH_VARARGS, "Create constant node"},
-    {"make_linear", py_make_linear, METH_VARARGS, "Create linear op node"},
-    {"make_log", py_make_log, METH_VARARGS, "Create log node"},
-    {"make_exp", py_make_exp, METH_VARARGS, "Create exp node"},
-    {"make_index", py_make_index, METH_VARARGS, "Create index node"},
-    {"make_add", py_make_add, METH_VARARGS, "Create add node"},
-    {"make_trace", py_make_trace, METH_VARARGS, "Create trace node"},
-    {"make_transpose", py_make_transpose, METH_VARARGS, "Create transpose node"},
-    {"make_hstack", py_make_hstack, METH_VARARGS,
-     "Create hstack node from list of expr capsules and n_vars (make_hstack([e1, "
-     "e2, ...], n_vars))"},
-    {"make_sum", py_make_sum, METH_VARARGS, "Create sum node"},
-    {"make_neg", py_make_neg, METH_VARARGS, "Create neg node"},
-    {"make_promote", py_make_promote, METH_VARARGS, "Create promote node"},
-    {"make_multiply", py_make_multiply, METH_VARARGS,
-     "Create elementwise multiply node"},
-    {"make_matmul", py_make_matmul, METH_VARARGS,
-     "Create matrix multiplication node (Z = X @ Y)"},
-    {"make_const_scalar_mult", py_make_const_scalar_mult, METH_VARARGS,
-     "Create constant scalar multiplication node (a * f(x))"},
-    {"make_const_vector_mult", py_make_const_vector_mult, METH_VARARGS,
-     "Create constant vector multiplication node (a ∘ f(x))"},
-    {"make_power", py_make_power, METH_VARARGS, "Create power node"},
-    {"make_prod", py_make_prod, METH_VARARGS, "Create prod node"},
-    {"make_prod_axis_zero", py_make_prod_axis_zero, METH_VARARGS,
-     "Create prod_axis_zero node"},
-    {"make_prod_axis_one", py_make_prod_axis_one, METH_VARARGS,
-     "Create prod_axis_one node"},
-    {"make_sin", py_make_sin, METH_VARARGS, "Create sin node"},
-    {"make_cos", py_make_cos, METH_VARARGS, "Create cos node"},
-    {"make_diag_vec", py_make_diag_vec, METH_VARARGS, "Create diag_vec node"},
-    {"make_tan", py_make_tan, METH_VARARGS, "Create tan node"},
-    {"make_sinh", py_make_sinh, METH_VARARGS, "Create sinh node"},
-    {"make_tanh", py_make_tanh, METH_VARARGS, "Create tanh node"},
-    {"make_asinh", py_make_asinh, METH_VARARGS, "Create asinh node"},
-    {"make_atanh", py_make_atanh, METH_VARARGS, "Create atanh node"},
-    {"make_broadcast", py_make_broadcast, METH_VARARGS, "Create broadcast node"},
-    {"make_entr", py_make_entr, METH_VARARGS, "Create entr node"},
-    {"make_logistic", py_make_logistic, METH_VARARGS, "Create logistic node"},
-    {"make_xexp", py_make_xexp, METH_VARARGS, "Create xexp node"},
-    {"make_left_matmul", py_make_left_matmul, METH_VARARGS,
-     "Create left matmul node (A @ f(x))"},
-    {"make_right_matmul", py_make_right_matmul, METH_VARARGS,
-     "Create right matmul node (f(x) @ A)"},
-    {"make_quad_form", py_make_quad_form, METH_VARARGS,
-     "Create quadratic form node (x' * Q * x)"},
-    {"make_quad_over_lin", py_make_quad_over_lin, METH_VARARGS,
-     "Create quad_over_lin node (sum(x^2) / y)"},
-    {"make_rel_entr", py_make_rel_entr, METH_VARARGS,
-     "Create rel_entr node: x * log(x/y) elementwise"},
-    {"make_rel_entr_vector_scalar", py_make_rel_entr_vector_scalar, METH_VARARGS,
-     "Create rel_entr node with vector first arg, scalar second arg"},
-    {"make_rel_entr_scalar_vector", py_make_rel_entr_scalar_vector, METH_VARARGS,
-     "Create rel_entr node with scalar first arg, vector second arg"},
-    {"get_expr_dimensions", py_get_expr_dimensions, METH_VARARGS,
-     "Get the dimensions (d1, d2) of an expression"},
-    {"get_expr_size", py_get_expr_size, METH_VARARGS,
-     "Get the total size of an expression"},
-    {"make_reshape", py_make_reshape, METH_VARARGS, "Create reshape atom"},
-    {"make_problem", py_make_problem, METH_VARARGS,
-     "Create problem from objective and constraints"},
-    {"problem_init_derivatives", py_problem_init_derivatives, METH_VARARGS,
-     "Initialize derivative structures"},
-    {"problem_objective_forward", py_problem_objective_forward, METH_VARARGS,
-     "Evaluate objective only"},
-    {"problem_constraint_forward", py_problem_constraint_forward, METH_VARARGS,
-     "Evaluate constraints only"},
-    {"problem_gradient", py_problem_gradient, METH_VARARGS,
-     "Compute objective gradient"},
-    {"problem_jacobian", py_problem_jacobian, METH_VARARGS,
-     "Compute constraint jacobian"},
-    {"get_jacobian", py_get_jacobian, METH_VARARGS,
-     "Get constraint jacobian without recomputing"},
-    {"problem_hessian", py_problem_hessian, METH_VARARGS,
-     "Compute Lagrangian Hessian"},
-    {"get_hessian", py_get_hessian, METH_VARARGS,
-     "Get Lagrangian Hessian without recomputing"},
-    {NULL, NULL, 0, NULL}};
-
-static struct PyModuleDef dnlp_module = {PyModuleDef_HEAD_INIT, "_diffengine", NULL,
-                                         -1, DNLPMethods};
-
-PyMODINIT_FUNC PyInit__diffengine(void)
-{
-    if (ensure_numpy() < 0) return NULL;
-    return PyModule_Create(&dnlp_module);
-}
diff --git a/python/problem/common.h b/python/problem/common.h
deleted file mode 100644
index 43e44c7..0000000
--- a/python/problem/common.h
+++ /dev/null
@@ -1,24 +0,0 @@
-#ifndef PROBLEM_COMMON_H
-#define PROBLEM_COMMON_H
-
-#define PY_SSIZE_T_CLEAN
-#include <Python.h>
-#include <numpy/arrayobject.h>
-
-#include "problem.h"
-
-/* Also need expr types for capsule handling */
-#include "../atoms/common.h"
-
-#define PROBLEM_CAPSULE_NAME "DNLP_PROBLEM"
-
-static void problem_capsule_destructor(PyObject *capsule)
-{
-    problem *prob = (problem *) PyCapsule_GetPointer(capsule, PROBLEM_CAPSULE_NAME);
-    if (prob)
-    {
-        free_problem(prob);
-    }
-}
-
-#endif /* PROBLEM_COMMON_H */
diff --git a/python/problem/constraint_forward.h b/python/problem/constraint_forward.h
deleted file mode 100644
index 9449db2..0000000
--- a/python/problem/constraint_forward.h
+++ /dev/null
@@ -1,54 +0,0 @@
-#ifndef PROBLEM_CONSTRAINT_FORWARD_H
-#define PROBLEM_CONSTRAINT_FORWARD_H
-
-#include "common.h"
-
-static PyObject *py_problem_constraint_forward(PyObject *self, PyObject *args)
-{
-    PyObject *prob_capsule;
-    PyObject *u_obj;
-    if (!PyArg_ParseTuple(args, "OO", &prob_capsule, &u_obj))
-    {
-        return NULL;
-    }
-
-    problem *prob =
-        (problem *) PyCapsule_GetPointer(prob_capsule, PROBLEM_CAPSULE_NAME);
-    if (!prob)
-    {
-        PyErr_SetString(PyExc_ValueError, "invalid problem capsule");
-        return NULL;
-    }
-
-    PyArrayObject *u_array =
-        (PyArrayObject *) PyArray_FROM_OTF(u_obj, NPY_DOUBLE, NPY_ARRAY_IN_ARRAY);
-    if (!u_array)
-    {
-        return NULL;
-    }
-
-    problem_constraint_forward(prob, (const double *) PyArray_DATA(u_array));
-    Py_DECREF(u_array);
-
-    PyObject *out = NULL;
-    if (prob->total_constraint_size > 0)
-    {
-        npy_intp size = prob->total_constraint_size;
-        out = PyArray_SimpleNew(1, &size, NPY_DOUBLE);
-        if (!out)
-        {
-            return NULL;
-        }
-        memcpy(PyArray_DATA((PyArrayObject *) out), prob->constraint_values,
-               size * sizeof(double));
-    }
-    else
-    {
-        npy_intp size = 0;
-        out = PyArray_SimpleNew(1, &size, NPY_DOUBLE);
-    }
-
-    return out;
-}
-
-#endif /* PROBLEM_CONSTRAINT_FORWARD_H */
diff --git a/python/problem/gradient.h b/python/problem/gradient.h
deleted file mode 100644
index 34b9b34..0000000
--- a/python/problem/gradient.h
+++ /dev/null
@@ -1,36 +0,0 @@
-#ifndef PROBLEM_GRADIENT_H
-#define PROBLEM_GRADIENT_H
-
-#include "common.h"
-
-static PyObject *py_problem_gradient(PyObject *self, PyObject *args)
-{
-    PyObject *prob_capsule;
-    if (!PyArg_ParseTuple(args, "O", &prob_capsule))
-    {
-        return NULL;
-    }
-
-    problem *prob =
-        (problem *) PyCapsule_GetPointer(prob_capsule, PROBLEM_CAPSULE_NAME);
-    if (!prob)
-    {
-        PyErr_SetString(PyExc_ValueError, "invalid problem capsule");
-        return NULL;
-    }
-
-    problem_gradient(prob);
-
-    npy_intp size = prob->n_vars;
-    PyObject *out = PyArray_SimpleNew(1, &size, NPY_DOUBLE);
-    if (!out)
-    {
-        return NULL;
-    }
-    memcpy(PyArray_DATA((PyArrayObject *) out), prob->gradient_values,
-           size * sizeof(double));
-
-    return out;
-}
-
-#endif /* PROBLEM_GRADIENT_H */
diff --git a/python/problem/hessian.h b/python/problem/hessian.h
deleted file mode 100644
index ea7f226..0000000
--- a/python/problem/hessian.h
+++ /dev/null
@@ -1,123 +0,0 @@
-#ifndef PROBLEM_HESSIAN_H
-#define PROBLEM_HESSIAN_H
-
-#include "common.h"
-
-/*
- * py_problem_hessian: Compute Lagrangian Hessian
- *
- * Args:
- *   prob_capsule: PyCapsule containing problem pointer
- *   obj_factor: Scaling factor for objective Hessian (double)
- *   lagrange: Array of Lagrange multipliers (numpy array, length =
- * total_constraint_size)
- *
- * Returns:
- *   Tuple of (data, indices, indptr, (m, n)) for scipy.sparse.csr_matrix
- */
-static PyObject *py_problem_hessian(PyObject *self, PyObject *args)
-{
-    PyObject *prob_capsule;
-    double obj_factor;
-    PyObject *lagrange_obj;
-
-    if (!PyArg_ParseTuple(args, "OdO", &prob_capsule, &obj_factor, &lagrange_obj))
-    {
-        return NULL;
-    }
-
-    problem *prob =
-        (problem *) PyCapsule_GetPointer(prob_capsule, PROBLEM_CAPSULE_NAME);
-    if (!prob)
-    {
-        PyErr_SetString(PyExc_ValueError, "invalid problem capsule");
-        return NULL;
-    }
-
-    /* Convert lagrange to contiguous C array */
-    PyArrayObject *lagrange_arr = (PyArrayObject *) PyArray_FROM_OTF(
-        lagrange_obj, NPY_DOUBLE, NPY_ARRAY_IN_ARRAY);
-    if (!lagrange_arr)
-    {
-        return NULL;
-    }
-
-    double *lagrange = (double *) PyArray_DATA(lagrange_arr);
-
-    /* Compute Hessian */
-    problem_hessian(prob, obj_factor, lagrange);
-
-    Py_DECREF(lagrange_arr);
-
-    /* Extract CSR components and return as tuple */
-    CSR_Matrix *H = prob->lagrange_hessian;
-    npy_intp nnz = H->nnz;
-    npy_intp n_plus_1 = H->n + 1;
-
-    PyObject *data = PyArray_SimpleNew(1, &nnz, NPY_DOUBLE);
-    PyObject *indices = PyArray_SimpleNew(1, &nnz, NPY_INT32);
-    PyObject *indptr = PyArray_SimpleNew(1, &n_plus_1, NPY_INT32);
-
-    if (!data || !indices || !indptr)
-    {
-        Py_XDECREF(data);
-        Py_XDECREF(indices);
-        Py_XDECREF(indptr);
-        return NULL;
-    }
-
-    /* Copy CSR data using memcpy for efficiency */
-    memcpy(PyArray_DATA((PyArrayObject *) data), H->x, nnz * sizeof(double));
-    memcpy(PyArray_DATA((PyArrayObject *) indices), H->i, nnz * sizeof(int));
-    memcpy(PyArray_DATA((PyArrayObject *) indptr), H->p, n_plus_1 * sizeof(int));
-
-    return Py_BuildValue("(OOO(ii))", data, indices, indptr, H->m, H->n);
-}
-
-static PyObject *py_get_hessian(PyObject *self, PyObject *args)
-{
-    PyObject *prob_capsule;
-    if (!PyArg_ParseTuple(args, "O", &prob_capsule))
-    {
-        return NULL;
-    }
-
-    problem *prob =
-        (problem *) PyCapsule_GetPointer(prob_capsule, PROBLEM_CAPSULE_NAME);
-    if (!prob)
-    {
-        PyErr_SetString(PyExc_ValueError, "invalid problem capsule");
-        return NULL;
-    }
-
-    if (!prob->lagrange_hessian)
-    {
-        PyErr_SetString(PyExc_RuntimeError,
-                        "hessian not initialized - call problem_hessian first");
-        return NULL;
-    }
-
-    CSR_Matrix *H = prob->lagrange_hessian;
-    npy_intp nnz = H->nnz;
-    npy_intp n_plus_1 = H->n + 1;
-
-    PyObject *data = PyArray_SimpleNew(1, &nnz, NPY_DOUBLE);
-    PyObject *indices = PyArray_SimpleNew(1, &nnz, NPY_INT32);
-    PyObject *indptr = PyArray_SimpleNew(1, &n_plus_1, NPY_INT32);
-
-    if (!data || !indices || !indptr)
-    {
-        Py_XDECREF(data);
-        Py_XDECREF(indices);
-        Py_XDECREF(indptr);
-        return NULL;
-    }
-
-    memcpy(PyArray_DATA((PyArrayObject *) data), H->x, nnz * sizeof(double));
-    memcpy(PyArray_DATA((PyArrayObject *) indices), H->i, nnz * sizeof(int));
-    memcpy(PyArray_DATA((PyArrayObject *) indptr), H->p, n_plus_1 * sizeof(int));
-
-    return Py_BuildValue("(OOO(ii))", data, indices, indptr, H->m, H->n);
-}
-
-#endif /* PROBLEM_HESSIAN_H */
diff --git a/python/problem/init_derivatives.h b/python/problem/init_derivatives.h
deleted file mode 100644
index cc6298d..0000000
--- a/python/problem/init_derivatives.h
+++ /dev/null
@@ -1,27 +0,0 @@
-#ifndef PROBLEM_INIT_DERIVATIVES_H
-#define PROBLEM_INIT_DERIVATIVES_H
-
-#include "common.h"
-
-static PyObject *py_problem_init_derivatives(PyObject *self, PyObject *args)
-{
-    PyObject *prob_capsule;
-    if (!PyArg_ParseTuple(args, "O", &prob_capsule))
-    {
-        return NULL;
-    }
-
-    problem *prob =
-        (problem *) PyCapsule_GetPointer(prob_capsule, PROBLEM_CAPSULE_NAME);
-    if (!prob)
-    {
-        PyErr_SetString(PyExc_ValueError, "invalid problem capsule");
-        return NULL;
-    }
-
-    problem_init_derivatives(prob);
-
-    Py_RETURN_NONE;
-}
-
-#endif /* PROBLEM_INIT_DERIVATIVES_H */
diff --git a/python/problem/jacobian.h b/python/problem/jacobian.h
deleted file mode 100644
index 1ba0a42..0000000
--- a/python/problem/jacobian.h
+++ /dev/null
@@ -1,105 +0,0 @@
-#ifndef PROBLEM_JACOBIAN_H
-#define PROBLEM_JACOBIAN_H
-
-#include "common.h"
-
-static PyObject *py_problem_jacobian(PyObject *self, PyObject *args)
-{
-    PyObject *prob_capsule;
-    if (!PyArg_ParseTuple(args, "O", &prob_capsule))
-    {
-        return NULL;
-    }
-
-    problem *prob =
-        (problem *) PyCapsule_GetPointer(prob_capsule, PROBLEM_CAPSULE_NAME);
-    if (!prob)
-    {
-        PyErr_SetString(PyExc_ValueError, "invalid problem capsule");
-        return NULL;
-    }
-
-    if (prob->n_constraints == 0)
-    {
-        /* Return empty CSR components */
-        npy_intp zero = 0;
-        npy_intp one = 1;
-        PyObject *data = PyArray_SimpleNew(1, &zero, NPY_DOUBLE);
-        PyObject *indices = PyArray_SimpleNew(1, &zero, NPY_INT32);
-        PyObject *indptr = PyArray_SimpleNew(1, &one, NPY_INT32);
-        ((int *) PyArray_DATA((PyArrayObject *) indptr))[0] = 0;
-        return Py_BuildValue("(OOO(ii))", data, indices, indptr, 0, prob->n_vars);
-    }
-
-    problem_jacobian(prob);
-
-    CSR_Matrix *jac = prob->jacobian;
-    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);
-        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));
-
-    return Py_BuildValue("(OOO(ii))", data, indices, indptr, jac->m, jac->n);
-}
-
-static PyObject *py_get_jacobian(PyObject *self, PyObject *args)
-{
-    PyObject *prob_capsule;
-    if (!PyArg_ParseTuple(args, "O", &prob_capsule))
-    {
-        return NULL;
-    }
-
-    problem *prob =
-        (problem *) PyCapsule_GetPointer(prob_capsule, PROBLEM_CAPSULE_NAME);
-    if (!prob)
-    {
-        PyErr_SetString(PyExc_ValueError, "invalid problem capsule");
-        return NULL;
-    }
-
-    if (!prob->jacobian)
-    {
-        PyErr_SetString(PyExc_RuntimeError,
-                        "jacobian not initialized - call problem_jacobian first");
-        return NULL;
-    }
-
-    CSR_Matrix *jac = prob->jacobian;
-    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);
-        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));
-
-    return Py_BuildValue("(OOO(ii))", data, indices, indptr, jac->m, jac->n);
-}
-
-#endif /* PROBLEM_JACOBIAN_H */
diff --git a/python/problem/make_problem.h b/python/problem/make_problem.h
deleted file mode 100644
index 6fc42d2..0000000
--- a/python/problem/make_problem.h
+++ /dev/null
@@ -1,66 +0,0 @@
-#ifndef PROBLEM_MAKE_H
-#define PROBLEM_MAKE_H
-
-#include "common.h"
-
-static PyObject *py_make_problem(PyObject *self, PyObject *args)
-{
-    PyObject *obj_capsule;
-    PyObject *constraints_list;
-    int verbose = 1;
-    if (!PyArg_ParseTuple(args, "OO|p", &obj_capsule, &constraints_list, &verbose))
-    {
-        return NULL;
-    }
-
-    expr *objective = (expr *) PyCapsule_GetPointer(obj_capsule, EXPR_CAPSULE_NAME);
-    if (!objective)
-    {
-        PyErr_SetString(PyExc_ValueError, "invalid objective capsule");
-        return NULL;
-    }
-
-    if (!PyList_Check(constraints_list))
-    {
-        PyErr_SetString(PyExc_TypeError, "constraints must be a list");
-        return NULL;
-    }
-
-    Py_ssize_t n_constraints = PyList_Size(constraints_list);
-    expr **constraints = NULL;
-    if (n_constraints > 0)
-    {
-        constraints = (expr **) malloc(n_constraints * sizeof(expr *));
-        if (!constraints)
-        {
-            PyErr_NoMemory();
-            return NULL;
-        }
-        for (Py_ssize_t i = 0; i < n_constraints; i++)
-        {
-            PyObject *c_capsule = PyList_GetItem(constraints_list, i);
-            constraints[i] =
-                (expr *) PyCapsule_GetPointer(c_capsule, EXPR_CAPSULE_NAME);
-            if (!constraints[i])
-            {
-                free(constraints);
-                PyErr_SetString(PyExc_ValueError, "invalid constraint capsule");
-                return NULL;
-            }
-        }
-    }
-
-    problem *prob =
-        new_problem(objective, constraints, (int) n_constraints, verbose);
-    free(constraints);
-
-    if (!prob)
-    {
-        PyErr_SetString(PyExc_RuntimeError, "failed to create problem");
-        return NULL;
-    }
-
-    return PyCapsule_New(prob, PROBLEM_CAPSULE_NAME, problem_capsule_destructor);
-}
-
-#endif /* PROBLEM_MAKE_H */
diff --git a/python/problem/objective_forward.h b/python/problem/objective_forward.h
deleted file mode 100644
index 51b1060..0000000
--- a/python/problem/objective_forward.h
+++ /dev/null
@@ -1,37 +0,0 @@
-#ifndef PROBLEM_OBJECTIVE_FORWARD_H
-#define PROBLEM_OBJECTIVE_FORWARD_H
-
-#include "common.h"
-
-static PyObject *py_problem_objective_forward(PyObject *self, PyObject *args)
-{
-    PyObject *prob_capsule;
-    PyObject *u_obj;
-    if (!PyArg_ParseTuple(args, "OO", &prob_capsule, &u_obj))
-    {
-        return NULL;
-    }
-
-    problem *prob =
-        (problem *) PyCapsule_GetPointer(prob_capsule, PROBLEM_CAPSULE_NAME);
-    if (!prob)
-    {
-        PyErr_SetString(PyExc_ValueError, "invalid problem capsule");
-        return NULL;
-    }
-
-    PyArrayObject *u_array =
-        (PyArrayObject *) PyArray_FROM_OTF(u_obj, NPY_DOUBLE, NPY_ARRAY_IN_ARRAY);
-    if (!u_array)
-    {
-        return NULL;
-    }
-
-    double obj_val =
-        problem_objective_forward(prob, (const double *) PyArray_DATA(u_array));
-
-    Py_DECREF(u_array);
-    return Py_BuildValue("d", obj_val);
-}
-
-#endif /* PROBLEM_OBJECTIVE_FORWARD_H */