"""PEtab-math translation for the ``observableFormula`` layer (issue #407, ADR-0035/0036).
Two production directions, both over ``petab``'s ``sympy``-backed math grammar
(``petab.v2.math``) -- PEtab math is a *specified* grammar, so we translate via its parsed
tree rather than a hand-rolled string tokenizer (ADR-0033 warned precisely against the
string approach -- operator precedence, the ``^`` power operator, and the
``ln``/``log10``/``sqrt`` spellings are where it would silently go wrong):
* :func:`bngl_body_to_petab_math` -- a BNGL function body -> a PEtab math expression
(the exporter's opt-in inlining mode, which generates the round-trip oracle). The body is
BNGL, not PEtab math, so it is read with BioNetGen's grammar (:mod:`._bngl_math`, #908),
printed as PEtab math with explicit parentheses, and guarded by a numeric check of the
printed formula against BioNetGen's reading of the body (:func:`_assert_matches_bngl`).
* :func:`compile_petab_formula` -- a PEtab math expression -> a vectorized ``numpy``
callable (the **measurement-model observation layer**, ADR-0036): the formula is evaluated
*post-simulation* over the output trajectory + the PSet, never by editing a model file.
This is the direction SBML import and the retrofitted BNGL-expression path turn on; it
superseded ADR-0035's *synthesis into the model file* (the ``PEtab-math -> BNGL function
body`` printer the importer once used to inject a ``begin functions`` entry).
``petab``/``sympy`` is the **optional runtime extra** ``pybnf[petab]`` -- imported lazily,
only on these expression paths. The bare-name ``observableFormula`` common case never
reaches this module and stays dependency-free + simulator-free (ADR-0019); an expression
import with ``petab`` absent raises a clear "install ``pybnf[petab]``" error, not an
``ImportError``.
**MVP scope (ADR-0035/0036).** Arithmetic over existing model entities (BNGL parameters /
observables / functions; SBML species / parameters) -- the surface a measurement model needs
(Boehm's quotient of sums is the worked fixture). A free symbol that is not a model entity is
an error; a PEtab ``observableParameter*``/``noiseParameter*`` per-measurement placeholder is
the deferred frontier and raises pointing here.
"""
import re
from ..printing import PybnfError
# A PEtab per-measurement placeholder symbol (``observableParameter1_*`` /
# ``noiseParameter1_*``): substituted per measurement row for scale/offset or a per-point
# noise value. It has no PyBNF analogue (PyBNF noise is per-observable, and there is no
# per-measurement observable scale/offset), so it is the deferred frontier, not a model
# entity (ADR-0035 / ADR-0033).
_PLACEHOLDER = re.compile(r'(?:observable|noise)Parameter\d')
def _require_petab_math():
"""The lazily-imported ``sympify_petab`` PEtab-math parser, or a pointed error.
``petab``/``sympy`` is the optional ``pybnf[petab]`` extra (ADR-0035): only the
expression path imports it. A missing install surfaces as a ``PybnfError`` naming the
extra, never a bare ``ImportError`` from deep in the call stack. Serialization is owned
by our own printers (:func:`._bngl_math.to_petab` for a BNGL body,
:func:`_petab_printer_cls` for a rewritten PEtab formula), not ``petab_math_str`` -- see
:func:`_petab_printer_cls` for why.
"""
try:
from petab.v2.math import sympify_petab
except ImportError as e:
raise PybnfError(
"An expression observableFormula needs the PEtab math translator, which is "
"the optional 'petab' extra. Install it with `pip install pybnf[petab]` (or "
"`uv pip install pybnf[petab]`). The bare-name observableFormula common case "
"(a model entity referenced by name) needs no translator and stays "
"dependency-free (ADR-0035, #407).") from e
return sympify_petab
# ---------------------------------------------------------------------------
# The translator pair
# ---------------------------------------------------------------------------
[docs]
def bngl_body_to_petab_math(body, entities, *, function_name=None, model_file=None):
"""Translate a BNGL function ``body`` to a PEtab math expression string.
The exporter's inlining mode (ADR-0035): a fitted **function** column emits its body
as ``observableFormula`` instead of the bare name. The body is BNGL, and BNGL does not
read arithmetic the way PEtab math does (``-k^2`` is ``(-k)^2`` and ``a^b^c`` is
``(a^b)^c`` in BNGL, the opposite in PEtab), so it is parsed with BioNetGen's grammar and
printed as PEtab math with explicit parentheses (:mod:`._bngl_math`, #908). A ``g()``
reference to another global function or an observable becomes the bare symbol ``g``
(PEtab math has no user zero-argument functions). Every free symbol is validated against
the model namespace (parameters u observables u functions), and the printed formula is
checked numerically against BioNetGen's reading of the body
(:func:`_assert_matches_bngl`) -- a wrong observableFormula is worse than a refused one
(ADR-0035). ``function_name``/``model_file`` only name the function in error messages.
Raises ``PybnfError`` on a missing ``petab`` extra, a malformed body, a body using an
operator BioNetGen's simulators refuse, an unknown free symbol, or a formula that
disagrees with the body; ``NotImplementedError`` on a construct PEtab math cannot
express exactly (``rint``, ``time()``, ``mratio``, ``TFUN``, a call with arguments to a
local function, a negative literal raised to a power, an ``if()`` whose condition is not
a comparison) and on a per-measurement placeholder symbol.
"""
from . import _bngl_math
sympify_petab = _require_petab_math()
where = (f"BNGL function '{function_name}'" if function_name else 'BNGL function body')
if model_file:
where += f" in model '{model_file}'"
try:
tree = _bngl_math.parse(
body, callable_names=set(entities.observable_names) | set(entities.function_names))
except _bngl_math.BnglBodyError as e:
raise PybnfError(f"Could not read the {where}, {body!r}: {e}.") from e
except NotImplementedError as e:
raise NotImplementedError(
f"The {where}, {body!r}, cannot be inlined as a PEtab observableFormula: {e}. "
f"Export without inline_functions to reference the function by name instead "
f"(ADR-0035, #908).") from e
petab_math = _bngl_math.to_petab(tree)
expr = _parse(sympify_petab, petab_math,
source=f'observableFormula printed for the {where}, {body!r}, which is')
_validate_symbols(expr, entities)
_assert_matches_bngl(tree, expr, petab_math, body, where)
return petab_math
[docs]
def inline_constants(formula, constants):
"""Substitute fixed-parameter constants into a PEtab math ``observableFormula``.
A PEtab ``observableFormula`` may reference a fixed parameter that lives only in the
PEtab parameters table, not in the model file (Boehm's ``specC17 = 0.107`` -- a
species-activity constant absent from the SBML; ADR-0037). Such a symbol resolves
against neither the model namespace nor the simulation trajectory, so the measurement
layer cannot evaluate it. Because it is *fixed*, substituting its numeric value is
exact: the measurement model then references only model entities + estimated
parameters and the model file stays unedited (ADR-0036).
``constants`` is the ``{name: value}`` map of fixed PEtab parameters **not** present
as model entities (a fixed parameter that *is* a model entity stays a symbol -- it
resolves as a model constant). Substitution + serialization go through ``sympy`` (never
a string tokenizer -- ADR-0033), guarded by the same numeric round-trip self-check as
the exporter. If no constant is a free symbol of ``formula`` the original string is
returned **verbatim** (the bare-name / model-only common case never reaches ``sympy``,
so the demo round trip is byte-stable).
Raises ``PybnfError`` on a missing ``petab`` extra, an unparseable formula, or a
substitution that fails its round-trip self-check.
"""
if not constants:
return formula
sympify_petab = _require_petab_math()
expr = _parse(sympify_petab, formula, source='observableFormula')
import sympy as sp
# Match by NAME: petab's parser tags symbols with assumptions (real/positive), so a plain
# ``sp.Symbol(name)`` is a different object -- substitute the actual free-symbol objects.
by_name = {str(s): s for s in expr.free_symbols}
present = {by_name[n]: v for n, v in constants.items() if n in by_name}
if not present:
return formula # nothing to inline -> carry the formula verbatim
subbed = expr.subs({sym: sp.Float(v) for sym, v in present.items()})
petab_math = _petab_printer_cls()().doprint(subbed)
_assert_round_trips(sympify_petab, subbed, petab_math, formula)
return petab_math
[docs]
def substitute_placeholders(formula, substitutions):
"""Substitute per-measurement placeholder symbols into a PEtab math ``formula`` (ADR-0044).
The sibling of :func:`inline_constants`, for the constant-per-observable placeholder
reduction: a PEtab ``observableParameter${n}_${id}`` / ``noiseParameter${n}_${id}``
placeholder whose measurements-table value is **constant across an observable's rows** is
not per-measurement at all -- it is a single scalar for that observable. Substituting it
reduces the placeholder to the existing per-observable machinery (the ADR-0036 measurement
layer for an ``observableFormula``; a sigma source for a ``noiseFormula``).
``substitutions`` maps a placeholder name to its token: a **number** is substituted as a
constant (``sp.Float``); a **parameter id** is substituted as a free symbol
(``sp.Symbol``), which resolves from the PSet at eval time (a nuisance free parameter,
ADR-0034). Substitution + serialization go through ``sympy`` (never a string tokenizer --
ADR-0033), guarded by the same numeric round-trip self-check as the exporter. A formula
with no substituted placeholder (the bare-name / no-placeholder common case) is returned
**verbatim** (it never reaches ``sympy``). A placeholder not in ``substitutions`` is left
in place (the caller validates it downstream).
Raises ``PybnfError`` on a missing ``petab`` extra, an unparseable formula, or a
substitution that fails its round-trip self-check.
"""
if not substitutions:
return formula
sympify_petab = _require_petab_math()
expr = _parse(sympify_petab, formula, source='formula')
import sympy as sp
# Match by NAME (petab tags symbols with assumptions, so a plain sp.Symbol(name) is a
# different object -- substitute the actual free-symbol objects), mirroring inline_constants.
by_name = {str(s): s for s in expr.free_symbols}
present = {by_name[n]: tok for n, tok in substitutions.items() if n in by_name}
if not present:
return formula # nothing to substitute -> carry the formula verbatim
repl = {}
for sym, token in present.items():
try:
repl[sym] = sp.Float(float(token)) # a numeric placeholder value -> a constant
except (TypeError, ValueError):
# A parameter id -> a free symbol (resolves from the PSet at eval). Build it through
# the petab parser so it carries the SAME assumptions (real/positive) re-parsing
# gives it; a plain sp.Symbol(token) would be a distinct object and the round-trip
# self-check below would false-reject the (correct) substitution.
repl[sym] = _parse(sympify_petab, token, source='observableParameters token')
subbed = expr.subs(repl)
petab_math = _petab_printer_cls()().doprint(subbed)
_assert_round_trips(sympify_petab, subbed, petab_math, formula)
return petab_math
[docs]
def inline_derived_symbols(formula, derived, *, observable_id=None):
"""Inline the derived SBML entities a measurement ``observableFormula`` references down to the
species/parameters they are defined over (#465 and #795, ADR-0036).
An SBML ``assignmentRule variable="X"`` makes ``X`` an algebraic function of other model
entities, recomputed every step -- never a simulation-output column and value-less, so the
measurement layer cannot resolve ``X`` *as a symbol* (#464). But its RHS *is* resolvable, and
SBML authors define exactly these convenience observables (the D2D ``Epo_cells := Epo_EpoRi +
dEpoi``), so ``observable: Epo_cells, formula: Epo_cells`` should just work. This substitutes
each referenced rule variable with its rule's RHS -- **recursively**, so a rule defined over
another rule resolves through (SBML allows forward references; resolution is by dependency,
not document order) -- leaving a formula over only species/parameters that the existing
:func:`compile_petab_formula` + the chain-rule gradient (:meth:`MeasurementModel.\
prediction_sensitivity`) handle unchanged: the rule's species sensitivities flow in automatically.
The second construct is an ``initialAssignment`` on a parameter or compartment (#795). SBML
lets one supersede the declared ``value``/``size``, so the attribute is a placeholder and the
entity has no value of its own. Inlining is what makes such a model work rather than merely
fail politely: an observable over Bertozzi's ``beta_N``, defined as ``(R0_ * gamma_) / N_``,
resolves to an expression over the three parameters, and if the fit estimates one of them the
measurement tracks the fit, because :meth:`MeasurementModel.materialize` resolves a symbol
from the parameter set before the constant snapshot. An initial assignment is only inlinable
when every entity it reads holds still for the whole simulation, since the substituted
expression is read at a measurement time while the assignment fixed a value at t=0; the
scanner applies that gate and records a non-inlinable definition as ``None``.
``derived`` maps every derived entity id to a ``DerivedSymbol`` carrying its ``kind``, its RHS
as a PEtab-math infix string (the stdlib ``_sbml`` scanner's serialization) or ``None``, and
the clause naming why a ``None`` cannot be inlined. A formula that references **no** derived
entity is returned **verbatim** (the common case never reaches sympy, so it stays byte-stable).
The substitution + serialization go through ``sympy`` (never a string tokenizer -- ADR-0033),
guarded by the same numeric round-trip self-check as the exporter.
Raises ``PybnfError`` on a missing ``petab`` extra, an unparseable formula/RHS, a reference to
a definition the scanner could not carry (``None`` RHS), a cyclic dependency, or a
substitution that fails its round-trip self-check.
"""
if not derived:
return formula
sympify_petab = _require_petab_math()
expr = _parse(sympify_petab, formula, source='observableFormula')
referenced = {str(s) for s in expr.free_symbols} & set(derived)
if not referenced:
return formula # no derived entity in the formula -> carry it verbatim
where = f"Measurement model '{observable_id}': " if observable_id else 'Measurement model: '
resolved = {} # name -> its fully derived-free sympy expr (memoized)
def resolve(name, path):
if name in resolved:
return resolved[name]
symbol = derived.get(name)
noun = ('initial assignment' if getattr(symbol, 'kind', None) == 'initial_assignment'
else 'assignment rule')
if name in path:
chain = ' -> '.join(path[path.index(name):] + (name,))
raise PybnfError(
f"{where}the SBML {noun} for '{name}' is cyclic ({chain}), so it cannot "
f"be resolved to a closed-form observable. (ADR-0036, #465 and #795.)")
rhs_infix = symbol.expression if symbol is not None else None
if rhs_infix is None:
reason = getattr(symbol, 'refusal', None) or (
f'its defining {noun} uses a construct this reader does not translate')
raise PybnfError(
f"{where}the observableFormula references '{name}', which the SBML model file "
f"defines in terms of other model entities and gives no value of its own. PyBNF "
f"cannot inline that definition here because {reason}. Write the observable "
f"directly over the species and parameters it is computed from instead. "
f"(ADR-0036, #465 and #795.)")
rhs = _parse(sympify_petab, rhs_infix, source=f"{noun} for {name!r}")
# Substitute by the ACTUAL free-symbol objects of this RHS parse (petab tags symbols with
# assumptions, so a plain sp.Symbol(name) would be a distinct object subs/diff ignore --
# the same gotcha inline_constants / substitute_placeholders guard).
sub = {s: resolve(str(s), path + (name,))
for s in rhs.free_symbols if str(s) in derived}
resolved[name] = rhs.subs(sub) if sub else rhs
return resolved[name]
by_name = {str(s): s for s in expr.free_symbols}
repl = {by_name[name]: resolve(name, ()) for name in referenced}
subbed = expr.subs(repl)
petab_math = _petab_printer_cls()().doprint(subbed)
_assert_round_trips(sympify_petab, subbed, petab_math, formula)
return petab_math
[docs]
def compile_objective_expression(formula, free_params, *, data_columns=(), backend='numpy'):
"""Compile a bring-your-own objective ``expression`` to ``(callable, ordered_names)``
(ADR-0050; the JAX backend is ADR-0059 item 2; data columns are the data-binding follow-up).
The fourth direction of the translator: the user writes a closed-form negative
log-likelihood (or cost) as PEtab math over the *declared free parameters* --
``objective = expression`` + ``expression = 0.5*((1 - x1)^2 + 100*(x2 - x1^2)^2)`` -- with
no model file (the inline analytical objective the #425 epic asks for). The sympy backend
is identical to :func:`compile_petab_formula`; the only differences are the validation
namespace and the wording: every free symbol must be a **declared free parameter** (so the
expression binds *by name* to the PSet, ADR-0050 §4 -- the bind-by-name footgun fix the
named-target slice deferred), not a model entity, and an unknown symbol is reported in those
terms (naming it). PEtab math uses ``^`` for exponentiation (not ``**``); a ``**`` surfaces
as a clear parse error here.
``free_params`` is the set/iterable of declared free-parameter names. ``data_columns`` is the
set/iterable of experimental-data column names *also* allowed as symbols when the objective
binds measurements (``data = curve.exp``; the data-binding follow-up): a data-bound expression
is a **per-observation** NLL contribution over the parameters *and* the row's data columns
(``(y - vmax*x/(km+x))^2`` references the free parameters ``vmax``/``km`` and the data columns
``x``/``y``), which the model evaluates per row and sums (``ExpressionModel.execute``). With no
bound data ``data_columns`` is empty and the contract is unchanged.
Returns the callable and the **sorted** free-symbol names it expects positionally
(``callable(*[binding[name] for name in ordered_names])``, where a parameter name binds to a
scalar and a data-column name to that column's array); a declared parameter the expression does
not reference is simply absent from the list (a likelihood flat in that direction).
``backend`` selects the lambdify target: ``'numpy'`` (the default) feeds the gradient-free
``score``-column path (``de`` / ``am`` / ``dream``); ``'jax'`` produces a JAX-traceable callable
so ``job_type = hmc`` can ``jax.grad`` the user's expression (ADR-0059 item 2). The parse,
symbol validation, and ``ordered_names`` are backend-independent (one sympy expression, one
sorted free-symbol set), so the numpy and JAX callables bind their arguments **identically** --
the bind-by-name contract holds across both, exactly as the prior families' numpy/JAX logpdfs
agree by construction.
Raises ``PybnfError`` on a missing ``petab`` extra, an unparseable expression, or a free
symbol that is neither a declared free parameter nor a bound data column."""
sympify_petab = _require_petab_math()
expr = _parse(sympify_petab, formula, source='objective expression')
params, columns = set(free_params), set(data_columns)
allowed = params | columns
column_hint = (f" or a bound data column ({sorted(columns)})" if columns else "")
_check_symbols(
expr, allowed,
unknown=lambda name: (
f"The objective expression references '{name}', which is not a declared free "
f"parameter{column_hint}. An inline 'objective = expression' may reference only "
f"parameters you declare (uniform_var / normal_var / a 'parameter:' record)"
f"{' and the columns of its bound data files' if columns else ''}; '{name}' is "
f"neither (ADR-0050). Either declare it as a free parameter or fix the typo."),
detail=f"Declared free parameters: {sorted(params)}." +
(f" Bound data columns: {sorted(columns)}." if columns else ""))
import sympy as sp
names = sorted(str(s) for s in expr.free_symbols)
func = sp.lambdify([sp.Symbol(n) for n in names], expr, modules=backend)
return func, names
# ---------------------------------------------------------------------------
# Shared helpers (parse, symbol validation, the BNGL <-> PEtab surface gap)
# ---------------------------------------------------------------------------
def _parse(sympify_petab, text, *, source):
"""Parse ``text`` to a sympy tree via PEtab's grammar (no evaluation, so the written
structure is preserved), turning a grammar error into a pointed ``PybnfError``."""
try:
return sympify_petab(text, evaluate=False)
except (ValueError, TypeError) as e:
raise PybnfError(
f"Could not parse the {source} {text!r} as PEtab math: {e}") from e
def _assert_matches_bngl(tree, expr, petab_math, body, where):
"""Refuse to emit a PEtab formula whose value differs from BioNetGen's reading of the body.
The export-inline safety net (ADR-0035, "a wrong measurement model is worse than a
refused one"; #908). ``tree`` is the body parsed with BioNetGen's grammar and ``expr`` is
PEtab's own parse of the printed ``petab_math``. At a fixed sequence of pseudo-random
points, taken in turn from three kinds (three quarters positive, spread over six decades,
so comparisons take both branches and a logarithm or a square root of a shifted quantity
finds its domain; all positive, where a fractional power of a quotient of several
symbols is defined; and small integers of either sign, where a negative base has a real
power), the tree is evaluated with BioNetGen's semantics in plain floating point
(:func:`._bngl_math.evaluate`), ``expr`` is evaluated the way the measurement layer will
evaluate it (``lambdify`` to numpy, so both sides are IEEE doubles), and the two must
agree. The two evaluations share nothing but the parse, so a printer that drops a
parenthesis PEtab needs, misspells a function, or mistranslates a comparison is caught
here, and so is a formula that petab's own reader would evaluate wrongly. (The guard this
replaced compared PEtab's parse of the body with PEtab's parse of the output, so a body
PEtab had misread in the first place passed it.)
A point where BioNetGen's value is undefined or already lost to double precision (a
domain error, a division by zero, an overflow, an underflow) is skipped. A point where
BioNetGen's value is finite and PEtab's is not, or where they differ, refuses the
formula; so does a body that is undefined at every point, since nothing could then be
checked.
"""
import random
import numpy as np
import sympy as sp
from . import _bngl_math
by_name = {str(s): s for s in expr.free_symbols}
names = sorted(_bngl_math.symbols(tree))
petab_names = sorted(by_name)
def refuse(detail):
raise PybnfError(
f"Refusing to emit the observableFormula {petab_math!r} for the {where}, "
f"{body!r}: {detail}, so emitting it would silently change the measurement model "
f"(a wrong observableFormula is worse than a refused one, ADR-0035, #908). Please "
f"report it, and export without inline_functions meanwhile.")
def unverifiable(detail):
raise PybnfError(
f"Could not check the observableFormula {petab_math!r} printed for the {where}, "
f"{body!r}, against the body, because {detail}. An unchecked formula is not "
f"emitted (ADR-0035, #908); export without inline_functions to reference the "
f"function by name.")
extra = sorted(set(by_name) - set(names))
if extra:
refuse(f"it names {extra}, which the body does not")
# petab parses log10(x) and log2(x) unevaluated as sympy's two-argument log(x, 10), which
# the numpy printer writes as numpy.log(x, 10) (the 10 lands in numpy's `out` slot and
# the call fails), so give lambdify a log that takes the base. Bases 10 and 2 use numpy's
# own log10 and log2, which are exact where log(x)/log(10) is not (log10(0.1) is -1.0,
# not -0.9999999999999998), so a power or a comparison fed by one is not misjudged.
def _log(x, base=None):
if base is None:
return np.log(x)
return {10: np.log10, 2: np.log2}.get(base, lambda v: np.log(v) / np.log(base))(x)
try:
func = sp.lambdify([by_name[n] for n in petab_names], expr,
modules=[{'log': _log}, 'numpy'])
except Exception as e: # noqa: BLE001 -- every failure here is a refusal, never a pass
unverifiable(f"sympy could not compile it for numerical evaluation "
f"({type(e).__name__}: {e})")
# Mixed signs alone left too few defined points for a Hill-type inverse such as
# (k^n*x/(r-x))^(1/n), and none for (-(2))^x, the form the refusal of (-2)^x recommends,
# so both were refused as uncheckable; hence the all-positive and the integer points
# (review of #908).
rng = random.Random(908)
agreed = 0
tries = 512
for i in range(tries):
if i % 3 == 0:
point = {n: rng.choice((1, 1, 1, -1)) * 10 ** rng.uniform(-3, 3) for n in names}
elif i % 3 == 1:
point = {n: 10 ** rng.uniform(-3, 3) for n in names}
else:
point = {n: float(rng.choice((1, 1, 1, -1)) * rng.randint(1, 4)) for n in names}
want = _bngl_math.evaluate(tree, point)
if want is None:
continue
try:
# numpy scalars, as the measurement layer passes arrays: Python floats would raise
# OverflowError where numpy (and a simulator) overflows to inf.
with np.errstate(all='ignore'):
got = float(func(*[np.float64(point[n]) for n in petab_names]))
except (ArithmeticError, TypeError, ValueError):
got = float('nan')
if not np.isfinite(got) or abs(got - want) > 1e-9 * max(1.0, abs(want)):
refuse(f"at {point} BioNetGen's value of the body is {want!r}, but petab's reading "
f"of the formula, evaluated as PyBNF's measurement layer evaluates it, is "
f"{got!r}")
agreed += 1
if agreed >= 8:
return
unverifiable(f'the body is undefined (a domain error, a division by zero or an overflow) '
f'at all but {agreed} of {tries} sample points, and 8 are needed')
def _assert_round_trips(sympify_petab, expr, petab_math, body):
"""Refuse to emit a PEtab serialization that does not parse back to ``expr``.
The safety net of the PEtab-to-PEtab rewrites (:func:`inline_constants`,
:func:`substitute_placeholders`, :func:`inline_derived_symbols`), whose input is already
PEtab math (ADR-0035, "a wrong measurement model is worse than a refused one"): re-parse
the emitted ``petab_math`` and assert it denotes the same function as ``expr``.
:func:`_petab_printer_cls` already parenthesizes the one petab serializer defect we know
of (the unparenthesized ``x ^ 1/2`` from a ``sqrt``); this guard is the standing tripwire
for *any* future serializer surprise, so corruption is always loud, never silent. It
cannot catch a misreading of the input, since both sides are read with PEtab's grammar;
that is why the BNGL export direction uses :func:`_assert_matches_bngl` instead (#908).
Equality is by **numeric sampling at several distinct positive points**, not symbolic
``simplify``/``equals``: petab floatifies literals (``sqrt`` parses back with a ``1.0/2.0``
Float exponent, not an exact ``Rational(1/2)``), and sympy treats Float-vs-exact powers
as undecidable -- so a symbolic test false-rejects the *correct* output. Positive points
keep ``sqrt``/``log`` real; multiple points rule out coincidental agreement (the corrupt
``z/2`` and ``sqrt(z)`` collide only at ``z=4``).
"""
if not _same_function(sympify_petab, expr, sympify_petab(petab_math, evaluate=False)):
raise PybnfError(
f"Refusing to emit the observableFormula {petab_math!r} for "
f"{body!r}: it does not parse back to the same function, so emitting it "
f"would silently corrupt the measurement model (a wrong observableFormula is "
f"worse than a refused one, ADR-0035). This indicates a PEtab math-serializer "
f"defect; please report it.")
def _same_function(sympify_petab, expr, other):
"""True iff ``expr`` and ``other`` evaluate equal at several distinct positive points."""
import sympy as sp
syms = sorted(expr.free_symbols | other.free_symbols, key=str)
agreed = 0
for k in range(1, 16):
subs = {s: sp.Rational(3 + 2 * k + 5 * i, 7) for i, s in enumerate(syms)}
try:
a, b = sp.N(expr.subs(subs)), sp.N(other.subs(subs))
except (TypeError, ValueError, ZeroDivisionError):
continue
if not (a.is_real and b.is_real and a.is_finite and b.is_finite):
continue # a domain/pole artifact at this point, not a translation error
if abs(float(a) - float(b)) > 1e-7 * max(1.0, abs(float(b))):
return False
agreed += 1
if agreed >= 4:
break
return True
def _namespace(entities):
"""The symbols an expression may reference: parameters u observables u functions.
Exactly the BNGL ``ParamList`` (ADR-0026): compartments, molecule types, and seed
species are not expression symbols, so a formula naming one is an error here.
"""
return (set(entities.parameters) | set(entities.observable_names)
| set(entities.function_names))
def _validate_symbols(expr, entities):
"""Assert every free symbol in ``expr`` is a known BNGL model entity.
An unknown symbol is an **error**, never a silent free parameter (ADR-0035). A PEtab
per-measurement placeholder (``observableParameter*`` / ``noiseParameter*``) is the
deferred frontier and raises ``NotImplementedError`` pointing at it; any other unknown
symbol raises ``PybnfError`` naming it and the model's entity sets.
"""
_check_symbols(
expr, _namespace(entities),
unknown=lambda name: (
f"The observableFormula references '{name}', which is not a parameter, "
f"observable, or function of the model. A measurement-model expression may "
f"only reference existing model entities (ADR-0035); an unknown symbol is an "
f"error, not a new free parameter."),
detail=(f"Model entities: parameters={sorted(entities.parameters)}; "
f"observables={sorted(entities.observable_names)}; "
f"functions={sorted(entities.function_names)}."))
def _check_symbols(expr, allowed, *, unknown, detail):
"""The shared free-symbol validator behind every translator direction (ADR-0035/0036).
Asserts each free symbol of ``expr`` is in ``allowed``. ``unknown(name)`` supplies the
direction-specific ``PybnfError`` summary (the BNGL translator and the numpy compiler
word it differently); ``detail`` is the namespace listing. A per-measurement placeholder
(``observableParameter*`` / ``noiseParameter*``) always raises the shared deferred-
frontier ``NotImplementedError`` -- one home for the placeholder boundary.
"""
for name in sorted(str(s) for s in expr.free_symbols):
if name in allowed:
continue
if _PLACEHOLDER.match(name):
raise NotImplementedError(
f"The observableFormula references a PEtab per-measurement placeholder "
f"'{name}' (an observableParameter*/noiseParameter* scale/offset or "
f"per-point noise value substituted per measurement row). PyBNF noise is "
f"per-observable and has no per-measurement observable scale/offset, so "
f"placeholders have no analogue and are the deferred frontier "
f"(ADR-0035 / ADR-0033 / ADR-0036, #407). This chunk translates arithmetic "
f"over existing model entities only.")
raise PybnfError(unknown(name), detail)
# Cached printer class. Defined lazily (it subclasses petab's printer) so this module
# imports with petab/sympy absent -- only the PEtab-to-PEtab rewrites build it. (The BNGL
# export direction prints with :func:`._bngl_math.to_petab` instead, #908.)
_PETAB_PRINTER = None
def _petab_printer_cls():
"""Build (once) and return the sympy-tree -> PEtab-math printer class.
Subclasses ``petab``'s own ``PetabStrPrinter`` (so functions/operators serialize
exactly as petab's validator expects) and overrides **only** ``_print_Pow``: petab
0.8.x's printer leaves a non-integer ``Rational`` exponent unparenthesized, so a
``sqrt`` (exponent ``1/2``) serializes as the precedence-unsafe ``x ^ 1/2`` -- which
re-parses as ``(x^1)/2`` and silently corrupts the measurement model. We parenthesize a
non-integer rational exponent (``x ^ (1/2)``), which both petab parses correctly and its
validator accepts. We own this rather than reverse ``petab_math_str`` so the forward
direction is precedence-safe, and :func:`_assert_round_trips` stands behind it as a
belt-and-suspenders check.
"""
global _PETAB_PRINTER
if _PETAB_PRINTER is not None:
return _PETAB_PRINTER
from petab.v2.math import PetabStrPrinter
class _PetabPrinter(PetabStrPrinter):
def _print_Pow(self, expr, rational=False):
base, exp = expr.as_base_exp()
str_base = self._print(base)
str_exp = self._print(exp)
if not base.is_Atom:
str_base = f'({str_base})'
# petab leaves a Rational atom like 1/2 unparenthesized ('x ^ 1/2' parses as
# (x^1)/2); parenthesize a non-integer rational (or any non-atom) exponent.
if not exp.is_Atom or (exp.is_Rational and not exp.is_Integer):
str_exp = f'({str_exp})'
return f'{str_base} ^ {str_exp}'
_PETAB_PRINTER = _PetabPrinter
return _PETAB_PRINTER
