Try out and potentially make use of SymPy's cse_diff.py module

[moorepants]

Found here: https://github.com/sympy/sympy/blob/master/sympy/simplify/_cse_diff.py

It is a simple swap to use it and it is include in SymPy 1.14. I tried it and all the opty tests passed. But the examples that use the gait2d model failed to differentiate. I don't think the assumptions of the original symbols are taken into account because I get this error:

 /home/moorepants/miniconda/envs/opty-dev/lib/python3.14/site-packages/sympy/printing/codeprinter.py(610)_print_not_supported()
    608     def _print_not_supported(self, expr):
    609         if self._settings.get('strict', False):
--> 610             raise PrintMethodNotImplementedError(
    611                 f"Unsupported by {type(self)}: {type(expr)}" +
    612                 "\nSet the printer option 'strict' to False in order to generate partially printed code."

ipdb> expr
re(x102)
ipdb> 

showing that the derivative of the foot penetration force is being treated as complex.

The Jacobian in cse_diff.py may be faster than the one in opty.

[moorepants]

diff for it to work:

diff --git a/opty/direct_collocation.py b/opty/direct_collocation.py
index ee42ccc..82535b8 100644
--- a/opty/direct_collocation.py
+++ b/opty/direct_collocation.py
@@ -7,6 +7,7 @@ import logging
 import numpy as np
 import sympy as sm
 from sympy.physics import mechanics as me
+from sympy.simplify._cse_diff import _forward_jacobian_norm_in_cse_out
 import cyipopt
 plt = sm.external.import_module('matplotlib.pyplot',
                                 import_kwargs={'fromlist': ['']},
@@ -2434,8 +2435,11 @@ class ConstraintCollocator(object):
         discrete_eom_matrix = sm.ImmutableDenseMatrix(self.discrete_eom)
         wrt_matrix = sm.ImmutableDenseMatrix([list(wrt)])
         if self._backend == 'cython':
-            symbolic_partials = _forward_jacobian(discrete_eom_matrix,
-                                                  wrt_matrix.T)
+            #symbolic_partials = _forward_jacobian(discrete_eom_matrix,
+                                                  #wrt_matrix.T)
+            r, j, _ = _forward_jacobian_norm_in_cse_out(discrete_eom_matrix,
+                                                        wrt_matrix.T)
+            symbolic_partials = (r, j)
         elif self._backend == 'numpy':
             symbolic_partials = discrete_eom_matrix.jacobian(wrt_matrix.T)
 

[moorepants]

This method could benefit by passing order='none' to cse() for performance gains.