Use finite differences in differentiate2c for implicit functions

python/nmodl/ode.py

@@ -643,15 +643,27 @@ def differentiate2c(expression, dependent_var, vars, prev_expressions=None):
     # differentiate w.r.t. x
     diff = expr.diff(x).simplify()
 
+    # could be something generic like f'(x), in which case we use finite differences
+    if needs_finite_differences(diff):
+        diff = (
+            transform_expression(diff, discretize_derivative)
+            .subs({finite_difference_step_variable(x): 1e-3})
+            .evalf()
+        )
+
+    # the codegen method does not like undefined function calls, so we extract
+    # them here
+    custom_fcts = {str(f.func): str(f.func) for f in diff.atoms(sp.Function)}
+
     # try to simplify expression in terms of existing variables
     # ignore any exceptions here, since we already have a valid solution
     # so if this further simplification step fails the error is not fatal
     try:
         # if expression is equal to one of the supplied vars, replace with this var
         # can do a simple string comparison here since a var cannot be further simplified
-        diff_as_string = sp.ccode(diff)
+        diff_as_string = sp.ccode(diff, user_functions=custom_fcts)
         for v in sympy_vars:
-            if diff_as_string == sp.ccode(sympy_vars[v]):
+            if diff_as_string == sp.ccode(sympy_vars[v], user_functions=custom_fcts):
                 diff = sympy_vars[v]
 
         # or if equal to rhs of one of the supplied equations, replace with lhs
@@ -672,4 +684,4 @@ def differentiate2c(expression, dependent_var, vars, prev_expressions=None):
         pass
 
     # return result as C code in NEURON format
-    return sp.ccode(diff.evalf())
+    return sp.ccode(diff.evalf(), user_functions=custom_fcts)

test/unit/ode/test_ode.py

@@ -4,6 +4,7 @@
 # SPDX-License-Identifier: Apache-2.0
 
 from nmodl.ode import differentiate2c, integrate2c
+import numpy as np
 
 import sympy as sp
 
@@ -100,6 +101,30 @@ def test_differentiate2c():
         "g",
     )
 
+    result = differentiate2c(
+        "-f(x)",
+        "x",
+        {},
+    )
+    # instead of comparing the expression as a string, we convert the string
+    # back to an expression and insert various functions
+    for function in [sp.sin, sp.exp, sp.tanh]:
+        for value in np.linspace(-5, 5, 100):
+            np.testing.assert_allclose(
+                float(
+                    sp.sympify(result)
+                    .subs(sp.Function("f"), function)
+                    .subs({"x": value})
+                    .evalf()
+                ),
+                float(
+                    -sp.Derivative(function("x"))
+                    .as_finite_difference(1e-3)
+                    .subs({"x": value})
+                    .evalf()
+                ),
+            )
+
 
 def test_integrate2c():