Use sbmlmath package for sympification of SBML math constructs

python/sdist/amici/de_model.py

@@ -409,7 +409,8 @@ def states(self) -> list[State]:
 
     def _process_sbml_rate_of(self) -> None:
         """Substitute any SBML-rateOf constructs in the model equations"""
-        rate_of_func = sp.core.function.UndefinedFunction("rateOf")
+        from sbmlmath import rate_of as rate_of_func
+
         species_sym_to_xdot = dict(
             zip(self.sym("x"), self.sym("xdot"), strict=True)
         )
@@ -449,25 +450,21 @@ def get_rate(symbol: sp.Symbol):
 
         # replace rateOf-instances in x0 by xdot equation
         for i_state in range(len(self.eq("x0"))):
-            if rate_ofs := self._eqs["x0"][i_state].find(rate_of_func):
-                self._eqs["x0"][i_state] = self._eqs["x0"][i_state].subs(
-                    {
-                        rate_of: get_rate(rate_of.args[0])
-                        for rate_of in rate_ofs
-                    }
-                )
+            new, replacement = self._eqs["x0"][i_state].replace(
+                rate_of_func, get_rate, map=True
+            )
+            if replacement:
+                self._eqs["x0"][i_state] = new
 
         # replace rateOf-instances in w by xdot equation
         #  here we may need toposort, as xdot may depend on w
         made_substitutions = False
         for i_expr in range(len(self.eq("w"))):
-            if rate_ofs := self._eqs["w"][i_expr].find(rate_of_func):
-                self._eqs["w"][i_expr] = self._eqs["w"][i_expr].subs(
-                    {
-                        rate_of: get_rate(rate_of.args[0])
-                        for rate_of in rate_ofs
-                    }
-                )
+            new, replacement = self._eqs["w"][i_expr].replace(
+                rate_of_func, get_rate, map=True
+            )
+            if replacement:
+                self._eqs["w"][i_expr] = new
                 made_substitutions = True
 
         if made_substitutions:

python/sdist/amici/import_utils.py

@@ -540,6 +540,8 @@ def _parse_heaviside_trigger(trigger: sp.Expr) -> sp.Expr:
             ]
         )
 
+    # TODO: x XOR y = (A & ~B) | (~A & B)
+    # TODO: x == y
     if isinstance(trigger, sp.And):
         return sp.Mul(*[_parse_heaviside_trigger(arg) for arg in trigger.args])
 
@@ -598,6 +600,7 @@ def _check_unsupported_functions(
         sp.functions.factorial,
         sp.functions.ceiling,
         sp.functions.floor,
+        sp.functions.tan,
         sp.functions.sec,
         sp.functions.csc,
         sp.functions.cot,