Symbolic

symbolic

Bounded scalar IR derivatives and safe Julia expression interoperability.

Derivatives describe the declared physical coordinates at regular domain points. They do not include initial-state sensitivities, solve an ODE or establish global identifiability. Native evaluation reuses the compiler's restricted evaluator.

DerivativeValues

Bases: NamedTuple

Borrow four native matrices without copying arrays or detaching graphs.

JuliaExport(model, state_names, parameter_names, observable_names, code, rational) dataclass

Bounded generated expressions and explicit external-runtime bindings.

t is the known independent time input. Initial vectors are supplied by the consumer in state_names order. Complete declarations remain in the companion PortableIR document. rational describes only the supported expression operations, not an external-tool identifiability conclusion.

SymbolicDerivatives(model, state_names, parameter_names, observable_names, expressions, state_dependencies, block_axes, _evaluate) dataclass

Ordered dense derivative artifacts for a bounded scalar-expression subset.

Shapes are (states, states), (states, parameters), (observables, states) and (observables, parameters). Structural dependency flags are independent of fitted parameter values. block_axes only repeats existing checked compiler declarations; no second separability analysis is performed.

evaluate(t, y, **parameters)

Evaluate all four matrices on native NumPy/JAX/Torch scalar bindings.

Returns:

Type	Description
DerivativeValues	Native matrices in the declared state/parameter/observable order.

Raises:

Type	Description
ValueError	Parameter or output layout differs.
TypeError	States do not have a native real floating dtype.

Source code in src/op_system/symbolic.py

def evaluate(self, t: object, y: object, **parameters: object) -> DerivativeValues:
    """Evaluate all four matrices on native NumPy/JAX/Torch scalar bindings.

    Returns:
        Native matrices in the declared state/parameter/observable order.

    Raises:
        ValueError: Parameter or output layout differs.
        TypeError: States do not have a native real floating dtype.
    """
    if set(parameters) != set(self.parameter_names):
        msg = "Derivative parameter names must match the declared order/layout"
        raise ValueError(msg)
    xp = array_namespace(y)
    if not xp.isdtype(getattr(y, "dtype", None), "real floating"):
        msg = "Derivative states require a native real floating array"
        raise TypeError(msg)
    values = self._evaluate(t, y, **parameters)
    if len(values.shape) != 1:
        msg = "Derivative parameters must be scalar; batch through native vmap"
        raise ValueError(msg)
    states, params, observed = (
        len(self.state_names),
        len(self.parameter_names),
        len(self.observable_names),
    )
    shapes = (
        (states, states),
        (states, params),
        (observed, states),
        (observed, params),
    )
    result = []
    offset = 0
    for shape in shapes:
        size = shape[0] * shape[1]
        result.append(xp.reshape(values[offset : offset + size], shape))
        offset += size
    return DerivativeValues(*result)

jvp(t, y, tangent, **parameters)

Apply the bounded dense state Jacobian.

Returns:

Type	Description
object	Native state-direction action; storage is dense, not matrix-free.

Raises:

Type	Description
ValueError	The tangent shape differs from the state layout.

Source code in src/op_system/symbolic.py

def jvp(
    self, t: object, y: object, tangent: object, **parameters: object
) -> object:
    """Apply the bounded dense state Jacobian.

    Returns:
        Native state-direction action; storage is dense, not matrix-free.

    Raises:
        ValueError: The tangent shape differs from the state layout.
    """
    matrix: Any = self.evaluate(t, y, **parameters).state
    if getattr(tangent, "shape", None) != (len(self.state_names),):
        msg = "State tangent shape differs from derivative layout"
        raise ValueError(msg)
    if array_namespace(tangent) is not array_namespace(y) or getattr(
        tangent, "dtype", None
    ) != getattr(y, "dtype", None):
        msg = "Tangent namespace/precision differs from the state"
        raise ValueError(msg)
    return matrix @ tangent

vjp(t, y, cotangent, **parameters)

Apply the transpose state Jacobian, retaining native derivative graphs.

Returns:

Type	Description
object	Native transpose action in the original state coordinates.

Raises:

Type	Description
ValueError	The cotangent shape differs from the state layout.

Source code in src/op_system/symbolic.py

def vjp(
    self, t: object, y: object, cotangent: object, **parameters: object
) -> object:
    """Apply the transpose state Jacobian, retaining native derivative graphs.

    Returns:
        Native transpose action in the original state coordinates.

    Raises:
        ValueError: The cotangent shape differs from the state layout.
    """
    matrix: Any = self.evaluate(t, y, **parameters).state
    if getattr(cotangent, "shape", None) != (len(self.state_names),):
        msg = "State cotangent shape differs from derivative layout"
        raise ValueError(msg)
    xp = array_namespace(y)
    if array_namespace(cotangent) is not xp or getattr(
        cotangent, "dtype", None
    ) != getattr(y, "dtype", None):
        msg = "Cotangent namespace/precision differs from the state"
        raise ValueError(msg)
    return xp.matrix_transpose(matrix) @ cotangent

export_julia(compiled, *, max_states=128)

Emit a bounded RHS function with explicit u/p/t binding order.

Only compiler-validated scalar IR emits code; original names never become Julia identifiers or executable interpolation. Initial states remain explicit caller inputs. The structural identity and ordered bindings are retained in the separate PortableIR document; emission is not a rank/tool certificate.

Returns:

Type	Description
JuliaExport	Generated code, source identity, ordered bindings and rational-subset flag.

Source code in src/op_system/symbolic.py

def export_julia(compiled: CompiledRhs, *, max_states: int = 128) -> JuliaExport:
    """Emit a bounded RHS function with explicit u/p/t binding order.

    Only compiler-validated scalar IR emits code; original names never become
    Julia identifiers or executable interpolation. Initial states remain explicit
    caller inputs. The structural identity and ordered bindings are retained in
    the separate PortableIR document; emission is not a rank/tool certificate.

    Returns:
        Generated code, source identity, ordered bindings and rational-subset flag.
    """
    equations, observed, parameters = _scalar_sources(compiled, max_states=max_states)
    bindings = {
        **{name: f"u[{i + 1}]" for i, name in enumerate(compiled.state_names)},
        **{name: f"p[{i + 1}]" for i, name in enumerate(parameters)},
        "t": "t",
    }
    rhs_text = (
        "function op_system_rhs(u, p, t)\n    return ["
        + ", ".join(_julia(expr, bindings) for expr in equations)
        + "]\nend\n"
    )
    observed_text = (
        "function op_system_observations(u, p, t)\n    return ["
        + ", ".join(_julia(expr, bindings) for _, expr in observed)
        + "]\nend\n"
    )
    rational = all(
        not (
            isinstance(node, Apply)
            and (
                node.op.startswith("np.")
                or (
                    node.op == "pow"
                    and (value := _exponent(node.args[1])) is not None
                    and value != int(value)
                )
            )
        )
        for expression in (*equations, *(expr for _, expr in observed))
        for node in _walk_scalar(expression)
    )
    return JuliaExport(
        StructuralModelIdentity.from_compiled(compiled),
        compiled.state_names,
        parameters,
        tuple(name for name, _ in observed),
        rhs_text + observed_text,
        rational,
    )

julia_rhs(compiled, *, max_states=128)

Return generated Julia code; prefer export_julia to retain its bindings.

Returns:

Type	Description
str	RHS and observation functions in the supported scalar subset.

Source code in src/op_system/symbolic.py

def julia_rhs(compiled: CompiledRhs, *, max_states: int = 128) -> str:
    """Return generated Julia code; prefer export_julia to retain its bindings.

    Returns:
        RHS and observation functions in the supported scalar subset.
    """
    return export_julia(compiled, max_states=max_states).code

symbolic_derivatives(compiled, *, max_states=128, max_entries=65536)

Export full state/parameter/observable derivatives for the supported IR.

Size checks precede dense derivative construction. Unsupported operations fail explicitly; numerical or AD fallbacks remain a caller choice.

Returns:

Type	Description
SymbolicDerivatives	Immutable ordered expressions and a native evaluator.

Raises:

Type	Description
ValueError	A dense artifact exceeds its declared size limits.

Source code in src/op_system/symbolic.py

def symbolic_derivatives(
    compiled: CompiledRhs, *, max_states: int = 128, max_entries: int = 65_536
) -> SymbolicDerivatives:
    """Export full state/parameter/observable derivatives for the supported IR.

    Size checks precede dense derivative construction. Unsupported operations
    fail explicitly; numerical or AD fallbacks remain a caller choice.

    Returns:
        Immutable ordered expressions and a native evaluator.

    Raises:
        ValueError: A dense artifact exceeds its declared size limits.
    """
    equations, observed, parameters = _scalar_sources(compiled, max_states=max_states)
    states = compiled.state_names
    entries = (len(states) + len(observed)) * (len(states) + len(parameters))
    if type(max_entries) is not int or max_entries < 1 or entries > max_entries:
        msg = "Symbolic dense derivative entry limit exceeded or invalid"
        raise ValueError(msg)
    groups = (
        tuple(_differentiate(expr, name) for expr in equations for name in states),
        tuple(_differentiate(expr, name) for expr in equations for name in parameters),
        tuple(_differentiate(expr, name) for _, expr in observed for name in states),
        tuple(
            _differentiate(expr, name) for _, expr in observed for name in parameters
        ),
    )
    expressions = tuple(tuple(unparse_ir(expr) for expr in group) for group in groups)
    flat = tuple(expr for group in groups for expr in group)
    evaluator = _make_eval_fn(
        state_names=states,
        aliases={},
        equations=tuple(unparse_ir(expr) for expr in flat),
        equations_ir=flat,
    )
    dependencies = tuple(
        tuple(name in free_symbols(expr) for name in states) for expr in equations
    )
    return SymbolicDerivatives(
        StructuralModelIdentity.from_compiled(compiled),
        states,
        parameters,
        tuple(name for name, _ in observed),
        expressions,
        dependencies,
        tuple(info.name for info in compiled.block_axes),
        evaluator,
    )