fix(symbolic): prevent Python execution in expression parsing

docs/symbolic/index.md

@@ -21,7 +21,7 @@
 - **`OutputEquation`** -- symbolic expression for sensor output `y` (e.g., `"x + noise"`)
 - **`compile_to_ode()`** -- lambdifies symbolic equations into a callable `ODEFunction` compatible with `gds-continuous`
 - **`linearize()`** -- computes Jacobian matrices (A, B, C, D) at an operating point
-- **Safe expression parsing** -- uses `sympy.parsing.sympy_parser.parse_expr`, never `eval`
+- **Restricted expression grammar** -- constructs SymPy expressions directly from validated syntax
 
 ## When to Use It
 
@@ -62,10 +62,10 @@ gds-control (pip install gds-domains)
 Symbolic expressions (strings)
     |
     v
-parse_expr()  -->  SymPy Expr objects
+Validated AST  -->  SymPy Expr objects
     |
     v
-compile_to_ode()  -->  ODEFunction (lambdified, numpy-backed)
+compile_to_ode()  -->  ODEFunction (lambdified, math-backed)
     |                       |
     v                       v
 linearize()           gds-continuous ODEModel
@@ -74,7 +74,39 @@ linearize()           gds-continuous ODEModel
 LinearSystem(A, B, C, D)   -->  eigenvalue analysis, controllability, etc.
 ```
 
-All expression parsing uses `sympy.parsing.sympy_parser.parse_expr` with a restricted transformation set -- arbitrary code execution is not possible.
+## Expression grammar and security
+
+ODE compilation, linearization (including output equations), and Hamiltonian
+state dynamics and Lagrangian parsing share a restricted expression grammar:
+
+- Integer and finite floating-point literals, including scientific notation.
+- Declared state, input, and parameter names; Hamiltonian expressions also
+  support `t` and generated costate names.
+- Constants `pi` and `E`, unless a declared variable uses the same name.
+- Parentheses, `+`, `-`, `*`, `/`, `**`, and unary `+`/`-`.
+- One-argument functions `sin`, `cos`, `tan`, `asin`, `acos`, `atan`, `sinh`,
+  `cosh`, `tanh`, `exp`, `sqrt`, and `Abs`; `log` takes one argument or an
+  optional second argument for the base.
+
+Use explicit multiplication (`2*x`) and `**` for powers. Unknown symbols,
+SymPy constructors such as `Symbol(...)`, attribute access, indexing, Python
+statements, arbitrary calls, keyword arguments, and argument unpacking raise
+`SymbolicError` when the expressions are consumed. Model construction performs
+structural validation; it does not parse expressions. `terminal_cost` is currently
+stored but not consumed by Hamiltonian derivation.
+
+The parser validates the complete Python AST and constructs SymPy objects directly;
+it does not evaluate the source string. Symbol names must be ASCII Python
+identifiers, excluding keywords and names starting with `__`. Code generation uses
+`lambdify` with dummy arguments and internally constructed expressions.
+`lambdify` itself generates and executes Python code, so arbitrary externally
+supplied SymPy objects are outside this expression-string security boundary.
+
+Expressions are limited to 4,096 characters, 512 AST nodes, depth 64, and integer
+literals of at most 256 bits. Constant simplification is deferred during parsing.
+These checks do not bound all subsequent symbolic differentiation, simplification,
+or numerical work. Services processing hostile input should enforce time and
+memory limits in an isolated worker.
 
 ## Quick Start
 

packages/gds-domains/README.md

@@ -20,3 +20,12 @@ pip install gds-domains[symbolic]                # + SymPy
 pip install gds-domains[nashpy]                  # + Nash equilibrium
 pip install gds-domains[all]                     # everything
 ```
+
+## Symbolic expressions
+
+Symbolic compilation and analysis accept a restricted mathematical grammar:
+declared variables, numeric literals, arithmetic, and approved math functions.
+Unsupported syntax raises `SymbolicError` when an expression is consumed.
+See the [expression grammar and security limits](https://dynamicalsystemsgroup.github.io/gds-core/symbolic/#expression-grammar-and-security)
+for supported functions, symbol naming rules, and resource limits. Parsing limits
+do not bound all subsequent symbolic computations.

packages/gds-domains/gds_domains/symbolic/_expressions.py

@@ -0,0 +1,167 @@
+"""Translate a bounded mathematical grammar to SymPy without evaluating Python."""
+
+from __future__ import annotations
+
+import ast
+import keyword
+import math
+from typing import Any
+
+from gds_domains.symbolic._compat import require_sympy
+from gds_domains.symbolic.errors import SymbolicError
+
+MAX_LENGTH = 4096
+MAX_NODES = 512
+MAX_DEPTH = 64
+MAX_INTEGER_BITS = 256
+
+# Explicit arities; no dynamic lookup of names on the SymPy module.
+_ARITIES = {
+    "sin": (1,),
+    "cos": (1,),
+    "tan": (1,),
+    "asin": (1,),
+    "acos": (1,),
+    "atan": (1,),
+    "sinh": (1,),
+    "cosh": (1,),
+    "tanh": (1,),
+    "exp": (1,),
+    "log": (1, 2),
+    "sqrt": (1,),
+    "Abs": (1,),
+}
+
+
+def make_symbols(names: list[str]) -> dict[str, Any]:
+    """Create symbols whose names can safely pass through code generation."""
+    require_sympy()
+    import sympy
+
+    for name in names:
+        if (
+            not name.isascii()
+            or not name.isidentifier()
+            or keyword.iskeyword(name)
+            or name.startswith("__")
+        ):
+            raise SymbolicError(
+                "Symbol names must be ASCII Python identifiers without "
+                "keywords or a leading double underscore"
+            )
+    return {name: sympy.Symbol(name) for name in names}
+
+
+def parse_expression(source: str, symbols: dict[str, Any]) -> Any:
+    """Parse only numbers, declared names, arithmetic and approved math calls.
+
+    Validate the entire tree before constructing any SymPy objects. Constructors
+    use evaluate=False to defer potentially expensive constant simplification.
+    These bounds do not guarantee limits on subsequent symbolic computations.
+    """
+    if len(source) > MAX_LENGTH:
+        raise SymbolicError(f"Expression exceeds {MAX_LENGTH} characters")
+    try:
+        tree = ast.parse(source.strip(), mode="eval")
+    except (SyntaxError, ValueError, RecursionError) as exc:
+        raise SymbolicError("Invalid mathematical expression syntax") from exc
+    if sum(1 for _ in ast.walk(tree)) > MAX_NODES:
+        raise SymbolicError(f"Expression exceeds {MAX_NODES} syntax nodes")
+
+    def validate(node: ast.AST, depth: int = 0) -> None:
+        if depth > MAX_DEPTH:
+            raise SymbolicError(f"Expression exceeds nesting depth {MAX_DEPTH}")
+        if isinstance(node, ast.Constant):
+            value = node.value
+            if type(value) is int:
+                if value.bit_length() > MAX_INTEGER_BITS:
+                    raise SymbolicError("Integer literal exceeds 256 bits")
+            elif type(value) is float:
+                if not math.isfinite(value):
+                    raise SymbolicError("Numeric literals must be finite")
+            else:
+                raise SymbolicError("Only real numeric literals are supported")
+        elif isinstance(node, ast.Name):
+            if node.id not in symbols and node.id not in {"pi", "E"}:
+                raise SymbolicError(f"Undeclared symbol: {node.id}")
+        elif isinstance(node, ast.UnaryOp) and isinstance(
+            node.op, (ast.UAdd, ast.USub)
+        ):
+            validate(node.operand, depth + 1)
+        elif isinstance(node, ast.BinOp) and isinstance(
+            node.op, (ast.Add, ast.Sub, ast.Mult, ast.Div, ast.Pow)
+        ):
+            validate(node.left, depth + 1)
+            validate(node.right, depth + 1)
+        elif isinstance(node, ast.Call) and isinstance(node.func, ast.Name):
+            if node.func.id not in _ARITIES:
+                raise SymbolicError(f"Unsupported math function: {node.func.id}")
+            if node.keywords or len(node.args) not in _ARITIES[node.func.id]:
+                raise SymbolicError(f"Invalid arguments for {node.func.id}")
+            for arg in node.args:
+                validate(arg, depth + 1)
+        else:
+            raise SymbolicError(f"Unsupported expression syntax: {type(node).__name__}")
+
+    validate(tree.body)
+    require_sympy()
+    import sympy
+
+    functions = {
+        "sin": sympy.sin,
+        "cos": sympy.cos,
+        "tan": sympy.tan,
+        "asin": sympy.asin,
+        "acos": sympy.acos,
+        "atan": sympy.atan,
+        "sinh": sympy.sinh,
+        "cosh": sympy.cosh,
+        "tanh": sympy.tanh,
+        "exp": sympy.exp,
+        "log": sympy.log,
+        "sqrt": sympy.sqrt,
+        "Abs": sympy.Abs,
+    }
+    constants = {"pi": sympy.pi, "E": sympy.E}
+
+    def build(node: ast.AST) -> Any:
+        if isinstance(node, ast.Constant):
+            return (
+                sympy.Integer(node.value)
+                if type(node.value) is int
+                else sympy.Float(node.value)
+            )
+        if isinstance(node, ast.Name):
+            return symbols[node.id] if node.id in symbols else constants[node.id]
+        if isinstance(node, ast.UnaryOp):
+            value = build(node.operand)
+            return (
+                value
+                if isinstance(node.op, ast.UAdd)
+                else sympy.Mul(-1, value, evaluate=False)
+            )
+        if isinstance(node, ast.BinOp):
+            left, right = build(node.left), build(node.right)
+            if isinstance(node.op, ast.Add):
+                return sympy.Add(left, right, evaluate=False)
+            if isinstance(node.op, ast.Sub):
+                return sympy.Add(
+                    left, sympy.Mul(-1, right, evaluate=False), evaluate=False
+                )
+            if isinstance(node.op, ast.Mult):
+                return sympy.Mul(left, right, evaluate=False)
+            if isinstance(node.op, ast.Div):
+                return sympy.Mul(
+                    left, sympy.Pow(right, -1, evaluate=False), evaluate=False
+                )
+            return sympy.Pow(left, right, evaluate=False)
+        if isinstance(node, ast.Call) and isinstance(node.func, ast.Name):
+            return functions[node.func.id](
+                *(build(arg) for arg in node.args), evaluate=False
+            )
+        raise AssertionError("Expression tree was not validated")
+
+    try:
+        return build(tree.body)
+    except (TypeError, ValueError, OverflowError, RecursionError) as exc:
+        raise SymbolicError("Unable to construct mathematical expression") from exc

packages/gds-domains/gds_domains/symbolic/compile.py

@@ -5,6 +5,7 @@
 from typing import TYPE_CHECKING, Any
 
 from gds_domains.symbolic._compat import require_sympy
+from gds_domains.symbolic._expressions import make_symbols, parse_expression
 
 if TYPE_CHECKING:
     from gds_continuous.types import ODEFunction
@@ -31,18 +32,17 @@ def compile_to_ode(
     input_names = [i.name for i in model.inputs]
 
     # Build symbol table
-    state_syms = {name: sympy.Symbol(name) for name in state_order}
-    input_syms = {name: sympy.Symbol(name) for name in input_names}
-    param_syms = {name: sympy.Symbol(name) for name in model.symbolic_params}
+    state_syms = make_symbols(state_order)
+    input_syms = make_symbols(input_names)
+    param_syms = make_symbols(model.symbolic_params)
 
     all_syms = {**state_syms, **input_syms, **param_syms}
 
-    # Parse expressions using safe parser (no eval, no builtins)
-    from sympy.parsing.sympy_parser import parse_expr
+    # Construct expressions from the restricted mathematical grammar.
 
     eq_map: dict[str, Any] = {}
     for eq in model.state_equations:
-        expr = parse_expr(eq.expr_str, local_dict=all_syms)
+        expr = parse_expression(eq.expr_str, all_syms)
         eq_map[eq.state_name] = expr
 
     # Build ordered RHS vector
@@ -60,7 +60,9 @@ def compile_to_ode(
         + [input_syms[n] for n in input_names]
         + [param_syms[n] for n in model.symbolic_params]
     )
-    rhs_lambda = sympy.lambdify(ordered_symbols, rhs_exprs, modules="math")
+    rhs_lambda = sympy.lambdify(
+        ordered_symbols, rhs_exprs, modules="math", dummify=True
+    )
 
     n_states = len(state_order)
 

packages/gds-domains/gds_domains/symbolic/elements.py

@@ -8,11 +8,11 @@
 class StateEquation(BaseModel, frozen=True):
     """Symbolic ODE right-hand side for a single state variable.
 
-    Declares: dx_i/dt = expr, where ``expr_str`` is a SymPy-parseable
-    string (e.g. ``"-k*x + u"``).
+    Declares: dx_i/dt = expr, where ``expr_str`` is a string in the
+    restricted mathematical grammar (e.g. ``"-k*x + u"``).
 
     The string form is R1-serializable. The sympy.Expr object is R3 —
-    reconstructed at lambdify time via ``sympy.sympify``.
+    constructed from validated syntax when the expression is consumed.
     """
 
     state_name: str

packages/gds-domains/gds_domains/symbolic/hamiltonian.py

@@ -17,6 +17,7 @@
 from pydantic import BaseModel
 
 from gds_domains.symbolic._compat import require_sympy
+from gds_domains.symbolic._expressions import make_symbols, parse_expression
 
 
 class HamiltonianSpec(BaseModel, frozen=True):
@@ -25,9 +26,9 @@ class HamiltonianSpec(BaseModel, frozen=True):
     Parameters
     ----------
     lagrangian
-        Running cost L(x, u, t) as a SymPy-parseable string.
+        Running cost L(x, u, t) using the restricted mathematical expression grammar.
     terminal_cost
-        Terminal cost Phi(x(T)) as a SymPy-parseable string.
+        Terminal cost Phi(x(T)) using the restricted mathematical expression grammar.
         Empty string means no terminal cost.
     control_bounds
         Lower and upper bounds for each control input.
@@ -100,14 +101,13 @@ def derive_hamiltonian(
     """
     require_sympy()
     import sympy
-    from sympy.parsing.sympy_parser import parse_expr
 
     # Build symbol tables
-    state_syms = {n: sympy.Symbol(n) for n in state_names}
+    state_syms = make_symbols(state_names)
     costate_names = [f"p_{n}" for n in state_names]
-    costate_syms = {n: sympy.Symbol(n) for n in costate_names}
-    input_syms = {n: sympy.Symbol(n) for n in input_names}
-    param_syms = {n: sympy.Symbol(n) for n in param_names}
+    costate_syms = make_symbols(costate_names)
+    input_syms = make_symbols(input_names)
+    param_syms = make_symbols(param_names)
     t_sym = sympy.Symbol("t")
 
     all_syms = {
@@ -122,12 +122,12 @@ def derive_hamiltonian(
     f_exprs = {}
     for name in state_names:
         if name in state_equations:
-            f_exprs[name] = parse_expr(state_equations[name], local_dict=all_syms)
+            f_exprs[name] = parse_expression(state_equations[name], all_syms)
         else:
             f_exprs[name] = sympy.Integer(0)
 
     # Parse Lagrangian L(x, u, t)
-    lagrangian = parse_expr(spec.lagrangian, local_dict=all_syms)
+    lagrangian = parse_expression(spec.lagrangian, all_syms)
 
     # Hamiltonian: H = L + p^T f
     hamiltonian = lagrangian
@@ -156,7 +156,9 @@ def derive_hamiltonian(
         + [input_syms[n] for n in input_names]
         + [param_syms[n] for n in param_names]
     )
-    rhs_lambda = sympy.lambdify(ordered_symbols, rhs_exprs, modules="math")
+    rhs_lambda = sympy.lambdify(
+        ordered_symbols, rhs_exprs, modules="math", dummify=True
+    )
 
     def augmented_ode(t: float, y: list[float], params: dict[str, Any]) -> list[float]:
         input_vals = [params.get(n, 0.0) for n in input_names]