Repository: PsiQ/bartiq

Issue:
Nested functions may fail to evaluate, depending on the order of functions in `custom_funcs` dict
## Describe the bug
Consider an expression like "f(g(x))". Substituting, in one go for `x`, `f` and `g` should result in a numerical value (as long as `f` and `g` produce numerical outputs for numerical outputs). However, this is not the case if the following two conditions are satisfied at the same time:

- In the input to `backend.substitute` method, `f` precedes `g` **and**
- Function assigned to `f` fails for nonnumerical inputs, in which case it is replaced by a dummy function of the same name.

## How to Reproduce
The following snippet reproduces the issue. Observe that `f` will fail for symbolic inputs.

```python
from bartiq.symbolics import sympy_backend

f = lambda x: int(x) + 1 # int(x) triggers a failure for symbolic expressions, which is needed here
g = lambda x: x ** 2

expr = sympy_backend.as_expression("f(g(x))")

result_1 = sympy_backend.substitute(expr, {"x": 3}, {"f": f, "g": g})
result_2 = sympy_backend.substitute(expr, {"x": 3}, {"g": g, "f": f})

print(result_1)  # Prints f(9), f failed to evaluate
print(result_2)  # Prints 10, this is correct
```

## Expected behavior
Both prints statements above produce the same output. The result of `backend.substitute` should not depend on the order of substitutions in input.

## Actual behavior
The `f` fails to evaluate in the first example above.

## Environment:**
This bug is OS/environment independent.

Discussion:

Patch:
diff --git a/src/bartiq/integrations/latex.py b/src/bartiq/integrations/latex.py
index b054f92..a865c98 100644
--- a/src/bartiq/integrations/latex.py
+++ b/src/bartiq/integrations/latex.py
@@ -18,7 +18,7 @@ from operator import attrgetter
 from qref.schema_v1 import ParamLinkV1, PortV1, ResourceV1, RoutineV1, SchemaV1
 from sympy import latex, symbols
 
-from ..symbolics.sympy_interpreter import parse_to_sympy
+from ..symbolics.sympy_backends import parse_to_sympy
 
 
 def routine_to_latex(routine: SchemaV1 | RoutineV1, show_non_root_resources: bool = True) -> str:
diff --git a/src/bartiq/symbolics/sympy_backends.py b/src/bartiq/symbolics/sympy_backends.py
index 2d82ea2..b32a1fd 100644
--- a/src/bartiq/symbolics/sympy_backends.py
+++ b/src/bartiq/symbolics/sympy_backends.py
@@ -95,6 +95,24 @@ def parse_to_sympy(expression: str, debug: bool = False) -> Expr:
     return parse(expression, interpreter=SympyInterpreter(debug=debug))
 
 
+def _sympify_function(func_name: str, func: Callable) -> type[sympy.Function]:
+    if not isinstance(func, sympy.Function):
+
+        def _eval_wrapper(cls, *args, **kwargs):
+            try:
+                return func(*args, **kwargs)
+            # The except clause here is intentionally broad, you never know what
+            # func can raise.
+            except Exception:
+                return None
+
+        sympy_func = type(func_name, (sympy.Function,), {"eval": classmethod(_eval_wrapper)})
+    else:
+        sympy_func = func
+
+    return sympy_func
+
+
 class SympyBackend:
 
     def __init__(self, parse_function: Callable[[str], Expr] = parse_to_sympy):
@@ -180,15 +198,11 @@ class SympyBackend:
                 f"Attempted to redefine the special function {func_name}; cannot define special functions."
             )
 
-        # Trying to evaluate a function which cannot be evaluated symbolically raises TypeError.
-        # This, however, is expected for certain functions (e.g. with conditions)
-        try:
-            return expr.replace(
-                lambda pattern: isinstance(pattern, SYMPY_USER_FUNCTION_TYPES) and str(type(pattern)) == func_name,
-                lambda match: function(*match.args),
-            )
-        except TypeError:
-            return expr
+        sympy_func = _sympify_function(func_name, function)
+        return expr.replace(
+            lambda pattern: isinstance(pattern, SYMPY_USER_FUNCTION_TYPES) and str(type(pattern)) == func_name,
+            lambda match: sympy_func(*match.args),
+        )
 
     def is_constant_int(self, expr: TExpr[Expr]):
         """Return True if a given expression represents a constant int and False otherwise."""


Test patch:
diff --git a/tests/symbolics/test_sympy_backend.py b/tests/symbolics/test_sympy_backend.py
index 105a40a..fe01774 100644
--- a/tests/symbolics/test_sympy_backend.py
+++ b/tests/symbolics/test_sympy_backend.py
@@ -147,3 +147,19 @@ def test_function_definition_succeeds_even_if_expression_becomes_constant(
     new_expr = backend.substitute(expr, variables, functions)
 
     assert backend.as_native(new_expr) == expected_native_result
+
+
+def test_defining_functions_is_invariant_under_input_order(backend):
+    def f(x):
+        return int(x) + 1
+
+    def g(x):
+        return x**2
+
+    expr = backend.as_expression("f(g(x))")
+
+    result_1 = backend.substitute(expr, {"x": 3}, {"f": f, "g": g})
+    result_2 = backend.substitute(expr, {"x": 3}, {"g": g, "f": f})
+
+    assert result_1 == 10
+    assert result_2 == 10
diff --git a/tests/symbolics/test_sympy_interpreter.py b/tests/symbolics/test_sympy_interpreter.py
index c334500..b49c596 100644
--- a/tests/symbolics/test_sympy_interpreter.py
+++ b/tests/symbolics/test_sympy_interpreter.py
@@ -386,6 +386,9 @@ PARSE_TEST_CASES = [
 
 @pytest.mark.parametrize("expression, expected_sympy_expression", PARSE_TEST_CASES)
 @pytest.mark.parametrize("parse", [parse_to_sympy, legacy_parse_to_sympy])
+@pytest.mark.filterwarnings("ignore:Legacy, pyparsing based sympy parser")
+@pytest.mark.filterwarnings(r"ignore:Using \^ operator to denote exponentiation is deprecated\.")
+@pytest.mark.filterwarnings(r"ignore:Results for using BigO with multiple #variables might be unreliable\.")
 def test_parse_to_sympy(expression, expected_sympy_expression, parse):
     """Tests for the sympy expression parser."""
     sympy_expression = parse(expression, debug=True)
@@ -395,3 +398,17 @@ def test_parse_to_sympy(expression, expected_sympy_expression, parse):
     new_expression = serialize_expression(sympy_expression)
     sympy_expression = parse(new_expression)
     assert sympy_expression == expected_sympy_expression
+
+
+def test_sympy_interpreter_warns_about_using_caret_sign_for_exponentiation():
+    expr = "x ^ 2"
+    with pytest.warns(match=r"Using \^ operator to denote exponentiation is deprecated\."):
+        _ = parse_to_sympy(expr)
+
+
+def test_legacy_parser_raises_deprecation_warning():
+    expr = "x + y"
+    with pytest.warns(match="Legacy, pyparsing based sympy parser"):
+        # Don't remove __wrapped__ here, we have to unwrap from lru_cache
+        # because otherwise we might miss the warning.
+        legacy_parse_to_sympy.__wrapped__(expr)
diff --git a/tests/test_analysis.py b/tests/test_analysis.py
index 219d6b6..02cf766 100644
--- a/tests/test_analysis.py
+++ b/tests/test_analysis.py
@@ -14,6 +14,7 @@
 
 
 import math
+import warnings
 
 import pytest
 import sympy
@@ -25,11 +26,6 @@ from bartiq.analysis import BigO, minimize
 @pytest.mark.parametrize(
     "expr,variable,expected",
     [
-        (
-            x**y + y**x + x * y + x**2 * y + 3 + x * y**2 + x + y + 1,
-            None,
-            BigO(x**y) + BigO(y**x) + BigO(x * y**2) + BigO(x**2 * y),
-        ),
         (
             x * y + x**2 * y + 3 + x * y**3 + x + y + 1,
             x,
@@ -45,14 +41,31 @@ from bartiq.analysis import BigO, minimize
             y,
             BigO(y),
         ),
+    ],
+)
+def test_BigO(expr, variable, expected):
+    assert BigO(expr, variable) == expected
+
+
+with warnings.catch_warnings():
+    warnings.filterwarnings("ignore", "Results for using BigO")
+    MULTIVARIATE_TEST_CASES = [
         (
             sympy.sympify("log(x) + y**2 + y"),
             None,
             BigO(sympy.sympify("log(x)")) + BigO(y**2),
         ),
-    ],
-)
-def test_BigO(expr, variable, expected):
+        (
+            x**y + y**x + x * y + x**2 * y + 3 + x * y**2 + x + y + 1,
+            None,
+            BigO(x**y) + BigO(y**x) + BigO(x * y**2) + BigO(x**2 * y),
+        ),
+    ]
+
+
+@pytest.mark.filterwarnings(r"ignore:Results for using BigO")
+@pytest.mark.parametrize("expr,variable,expected", MULTIVARIATE_TEST_CASES)
+def test_multivariate_BigO(expr, variable, expected):
     assert BigO(expr, variable) == expected
 
 
@@ -64,6 +77,7 @@ def test_BigO_throws_warning_for_multiple_variables():
         BigO(x**y + y**x + x * y + x**2 * y + 3 + x * y**2 + x + y + 1)
 
 
+@pytest.mark.filterwarnings(r"ignore:Results for using BigO")
 def test_adding_BigO_expressions():
     assert BigO(x) + BigO(x) == BigO(x)
     assert BigO(x) * BigO(x) == BigO(x**2)
@@ -158,7 +172,7 @@ def test_minimize_gradient_descent(
             {
                 "method": "L-BFGS-B",
                 "tol": 1e-6,
-                "options": {"disympy": False},
+                "options": {"disp": False},
             },
             math.pi,
             -1.0,
@@ -175,7 +189,7 @@ def test_minimize_gradient_descent(
             {
                 "method": "Nelder-Mead",
                 "tol": 1e-6,
-                "options": {"disympy": False},
+                "options": {"disp": False},
             },
             0.0,
             0.0,