Diffengine extractor: DNLP solving with differentiation engine

Squashed 163 commits from diffengine-extractor branch. Includes:
- DNLP (Disciplined Nonlinear Programming) framework
- NLP solver interfaces (IPOPT, COPT, Knitro, Uno)
- Differentiation engine with sparse Jacobian/Hessian support
- dnlp2smooth canonicalizers for all supported atoms
- NLP solving chain and problem reduction pipeline
- Comprehensive NLP test suite

Author: Transurgeon

CLAUDE.md

@@ -1,24 +1,36 @@
+# CLAUDE.md
+
+This file provides guidance to Claude Code (claude.ai/code) when working with code in this repository.
+
 # CVXPY Development Guide
 
 ## Quick Reference
 
 ### Commands
 ```bash
-# Install in development mode
+# Install in development mode (includes C++ extensions)
 pip install -e .
 
 # Install pre-commit hooks (required)
 pip install pre-commit && pre-commit install
 
+# Run linter manually
+ruff check cvxpy/ --fix
+
 # Run all tests
 pytest cvxpy/tests/
 
 # Run specific test
 pytest cvxpy/tests/test_atoms.py::TestAtoms::test_norm_inf
+
+# Run benchmarks (uses asv - airspeed velocity)
+cd benchmarks && pip install -e . && asv run
 ```
 
 ## Code Style
 
+- **Python >= 3.11** required
+- **Linter**: Ruff (enforced via pre-commit hooks; auto-fixes code but fails the commit if fixes were needed). Note: `__init__.py` files are excluded from linting.
 - **Line length**: 100 characters
 - **IMPORTANT: IMPORTS AT THE TOP** of files - circular imports are the only exception
 - **IMPORTANT:** Add Apache 2.0 license header to all new files
@@ -126,14 +138,13 @@ Check with `problem.is_dpp()`. See `cvxpy/utilities/scopes.py` for implementatio
 Location: `cvxpy/atoms/` or `cvxpy/atoms/elementwise/`
 
 ```python
-from typing import Tuple
 from cvxpy.atoms.atom import Atom
 
 class my_atom(Atom):
     def __init__(self, x) -> None:
         super().__init__(x)
 
-    def shape_from_args(self) -> Tuple[int, ...]:
+    def shape_from_args(self) -> tuple[int, ...]:
         return self.args[0].shape
 
     def sign_from_args(self) -> Tuple[bool, bool]:

cvxpy/atoms/__init__.py

@@ -163,6 +163,12 @@
     ptp
 ]
 
+NON_SMOOTH_ATOMS = [
+    abs,
+    maximum,
+    minimum,
+]
+
 # DGP atoms whose Dgp2Dcp canonicalization produces ExpCone-requiring DCP atoms.
 GP_EXP_ATOMS = [
     AddExpression,

cvxpy/atoms/affine/transpose.py

@@ -116,6 +116,7 @@ def graph_implementation(
         """
         return (lu.transpose(arg_objs[0], self.axes), [])
 
+
 def permute_dims(expr, axes: List[int]):
     """Permute the dimensions of the expression.
 

cvxpy/atoms/elementwise/entr.py

@@ -97,3 +97,6 @@ def _domain(self) -> List[Constraint]:
         """Returns constraints describing the domain of the node.
         """
         return [self.args[0] >= 0]
+    
+    def point_in_domain(self):
+        return np.ones(self.shape)

cvxpy/atoms/elementwise/kl_div.py

@@ -98,3 +98,6 @@ def _domain(self) -> List[Constraint]:
         """Returns constraints describing the domain of the node.
         """
         return [self.args[0] >= 0, self.args[1] >= 0]
+    
+    def point_in_domain(self, argument=0):
+        return np.ones(self.args[argument].shape)

cvxpy/atoms/elementwise/log.py

@@ -105,4 +105,8 @@ def _domain(self) -> List[Constraint]:
         """Returns constraints describing the domain of the node.
         """
         return [self.args[0] >= 0]
-    
+
+    def point_in_domain(self) -> np.ndarray:
+        """Returns a point in the domain of the node.
+        """
+        return np.ones(self.shape)

cvxpy/atoms/elementwise/rel_entr.py

@@ -99,4 +99,6 @@ def _domain(self):
         """Returns constraints describing the domain of the node.
         """
         return [self.args[0] >= 0, self.args[1] >= 0]
-    
+
+    def point_in_domain(self, argument=0):
+        return np.ones(self.args[argument].shape)

cvxpy/atoms/elementwise/xexp.py

@@ -96,4 +96,4 @@ def _grad(self, values):
     def _domain(self) -> List[Constraint]:
         """Returns constraints describing the domain of the node.
         """
-        return [self.args[0] >= 0]
+        return [self.args[0] >= 0]

cvxpy/reductions/chain.py

@@ -1,3 +1,5 @@
+import time
+
 from cvxpy import settings as s
 from cvxpy.reductions.reduction import Reduction
 
@@ -73,7 +75,11 @@ def apply(self, problem, verbose: bool = False):
         for r in self.reductions:
             if verbose:
                 s.LOGGER.info('Applying reduction %s', type(r).__name__)
+            start = time.perf_counter()
             problem, inv = r.apply(problem)
+            elapsed = time.perf_counter() - start
+            if verbose:
+                s.LOGGER.info('  %s took %.4f seconds', type(r).__name__, elapsed)
             inverse_data.append(inv)
         return problem, inverse_data
 

cvxpy/reductions/dcp2cone/cone_matrix_stuffing.py

@@ -322,6 +322,53 @@ def split_solution(self, sltn, active_vars=None):
                         value, var.shape, order='F')
         return sltn_dict
 
+    def apply_restruct_mat(self, restruct_mat, restruct_mat_op):
+        """Apply restructuring matrix to parametric tensor A.
+
+        Parameters
+        ----------
+        restruct_mat : list
+            List of sparse matrices or linear operators (unused for this path).
+        restruct_mat_op : LinearOperator
+            Block diagonal linear operator for restructuring.
+
+        Returns
+        -------
+        ParamConeProg
+            New program with restructured A tensor.
+        """
+        if restruct_mat_op is not None:
+            unspecified, _ = np.divmod(self.A.shape[0] * self.A.shape[1],
+                                       restruct_mat_op.shape[1], dtype=np.int64)
+            reshaped_A = self.A.reshape(restruct_mat_op.shape[1],
+                                        unspecified, order='F').tocsr()
+            restructured_A = restruct_mat_op(reshaped_A).tocoo()
+            # Because of a bug in scipy versions < 1.20, `reshape`
+            # can overflow if indices are int32s.
+            restructured_A.row = restructured_A.row.astype(np.int64)
+            restructured_A.col = restructured_A.col.astype(np.int64)
+            restructured_A = restructured_A.reshape(
+                np.int64(restruct_mat_op.shape[0]) * (np.int64(self.x.size) + 1),
+                self.A.shape[1], order='F')
+        else:
+            restructured_A = self.A
+        return ParamConeProg(
+            self.q,
+            self.x,
+            restructured_A,
+            self.variables,
+            self.var_id_to_col,
+            self.constraints,
+            self.parameters,
+            self.param_id_to_col,
+            P=self.P,
+            formatted=True,
+            lower_bounds=self.lower_bounds,
+            upper_bounds=self.upper_bounds,
+            lb_tensor=self.lb_tensor,
+            ub_tensor=self.ub_tensor,
+        )
+
     def split_adjoint(self, del_vars=None):
         """Adjoint of split_solution.
         """
@@ -339,6 +386,65 @@ def split_adjoint(self, del_vars=None):
         return var_vec
 
 
+def lower_and_order_constraints(constraints):
+    """Lower equality/inequality constraints and reorder by cone type.
+
+    Converts Equality -> Zero, Inequality/NonPos -> NonNeg, and normalizes
+    SOC/PowCone/ExpCone axes. Returns constraints ordered as:
+    Zero, NonNeg, SOC, PSD, ExpCone, PowCone3D, PowConeND.
+
+    Parameters
+    ----------
+    constraints : list
+        The problem constraints to lower and reorder.
+
+    Returns
+    -------
+    ordered_cons : list
+        The lowered and reordered constraints.
+    cons_id_map : dict
+        Mapping from constraint id to constraint id (identity map).
+    """
+    cons = []
+    for con in constraints:
+        if isinstance(con, Equality):
+            con = lower_equality(con)
+        elif isinstance(con, Inequality):
+            con = lower_ineq_to_nonneg(con)
+        elif isinstance(con, NonPos):
+            con = nonpos2nonneg(con)
+        elif isinstance(con, SOC) and con.axis == 1:
+            con = SOC(con.args[0], con.args[1].T, axis=0,
+                      constr_id=con.constr_id)
+        elif isinstance(con, PowCone3D) and con.args[0].ndim > 1:
+            x, y, z = con.args
+            alpha = con.alpha
+            con = PowCone3D(x.flatten(order='F'),
+                            y.flatten(order='F'),
+                            z.flatten(order='F'),
+                            alpha.flatten(order='F'),
+                            constr_id=con.constr_id)
+        elif isinstance(con, PowConeND) and con.axis == 1:
+            alpha = con.alpha.T
+            W = con.W.T
+            con = PowConeND(W, con.z.flatten(order='F'),
+                            alpha,
+                            axis=0,
+                            constr_id=con.constr_id)
+        elif isinstance(con, ExpCone) and con.args[0].ndim > 1:
+            x, y, z = con.args
+            con = ExpCone(x.flatten(order='F'), y.flatten(order='F'), z.flatten(order='F'),
+                          constr_id=con.constr_id)
+        cons.append(con)
+    # Reorder constraints to Zero, NonNeg, SOC, PSD, EXP, PowCone3D, PowConeND
+    constr_map = group_constraints(cons)
+    ordered_cons = constr_map[Zero] + constr_map[NonNeg] + \
+        constr_map[SOC] + constr_map[PSD] + constr_map[ExpCone] + \
+        constr_map[PowCone3D] + constr_map[PowConeND]
+    cons_id_map = {con.id: con.id for con in ordered_cons}
+    return ordered_cons, cons_id_map
+
+
 class ConeMatrixStuffing(MatrixStuffing):
     """Construct matrices for linear cone problems.
 
@@ -380,51 +486,20 @@ def stuffed_objective(self, problem, extractor):
 
     def apply(self, problem):
         inverse_data = InverseData(problem)
-        # Lower equality and inequality to Zero and NonNeg.
-        cons = []
-        for con in problem.constraints:
-            if isinstance(con, Equality):
-                con = lower_equality(con)
-            elif isinstance(con, Inequality):
-                con = lower_ineq_to_nonneg(con)
-            elif isinstance(con, SOC) and con.axis == 1:
-                con = SOC(con.args[0], con.args[1].T, axis=0,
-                          constr_id=con.constr_id)
-            elif isinstance(con, PowCone3D) and con.args[0].ndim > 1:
-                x, y, z = con.args
-                alpha = con.alpha
-                con = PowCone3D(x.flatten(order='F'),
-                                y.flatten(order='F'),
-                                z.flatten(order='F'),
-                                alpha.flatten(order='F'),
-                                constr_id=con.constr_id)
-            elif isinstance(con, PowConeND) and con.axis == 1:
-                alpha = con.alpha.T
-                W = con.W.T
-                con = PowConeND(W, con.z.flatten(order='F'),
-                                alpha,
-                                axis=0,
-                                constr_id=con.constr_id)
-            elif isinstance(con, ExpCone) and con.args[0].ndim > 1:
-                x, y, z = con.args
-                con = ExpCone(x.flatten(order='F'), y.flatten(order='F'), z.flatten(order='F'),
-                              constr_id=con.constr_id)
-            cons.append(con)
+        # Lower equality and inequality constraints and reorder by cone type.
+        ordered_cons, cons_id_map = lower_and_order_constraints(problem.constraints)
+        inverse_data.cons_id_map = cons_id_map
+        inverse_data.constraints = ordered_cons
+
         # Need to check that intended canonicalization backend still works.
-        lowered_con_problem = problem.copy([problem.objective, cons])
+        lowered_con_problem = problem.copy(
+            [problem.objective, ordered_cons])
         canon_backend = get_canon_backend(lowered_con_problem, self.canon_backend)
         # Form the constraints
         extractor = CoeffExtractor(inverse_data, canon_backend)
         params_to_P, params_to_c, flattened_variable = self.stuffed_objective(
             problem, extractor)
-        # Reorder constraints to Zero, NonNeg, SOC, PSD, EXP, PowCone3D, PowConeND
-        constr_map = group_constraints(cons)
-        ordered_cons = constr_map[Zero] + constr_map[NonNeg] + \
-            constr_map[SOC] + constr_map[PSD] + constr_map[ExpCone] + \
-            constr_map[PowCone3D] + constr_map[PowConeND]
-        inverse_data.cons_id_map = {con.id: con.id for con in ordered_cons}
 
-        inverse_data.constraints = ordered_cons
         # Batch expressions together, then split apart.
         expr_list = [arg for c in ordered_cons for arg in c.args]
         params_to_problem_data = extractor.affine(expr_list)