"""PEtab v2 exporter: a PyBNF/BNGL job -> PEtab v2 artifacts (#407/#423; ADR-0025/27/28).
The **exporter-first** direction of the PEtab interop (ADR-0025): a working PyBNF
BNGL job and a native ``.conf`` are *read* and serialized to a PEtab v2 problem
(``parameters.tsv`` / ``observables.tsv`` / ``measurements.tsv`` / ``conditions.tsv`` /
``experiments.tsv`` / ``problem.yaml`` + a PEtab-clean copy of the model), rather than
*generating* BNGL from a declarative PEtab spec (the harder importer direction, deferred).
The reverse asset mappings live beside their importer twins
(``parameters.petab_parameter_row``, ``observables.petab_observable_row``,
``measurements.measurement_rows_from_data``, ``conditions.build_*``); this module is the
*disposable* glue: it reads the job (the stdlib ``ploop`` config parser, a focused BNGL
block reader, and :class:`pybnf.data.Data` for the ``.exp``) and writes the files.
**Why a function is the measurement model.** A fitted ``.exp`` column matches a BNGL
**observable** *or* a **function** (PyBNF forces ``print_functions=>1``), and the
function is usually the measurement model. So an observable column exports to
``observableFormula = <name>`` and a function column to ``observableFormula =
<name>`` too -- always the *bare model name*, with the function carried verbatim in
the model file (ADR-0025). PEtab ids are prefixed (``obs_``/``func_`` for
observables, the unprefixed model name for parameters) to keep the PEtab-id namespace
disjoint from the model-entity namespace.
**New-era only (ADR-0028/0031): export is transcription.** PEtab v2 interop is a new-era
feature, and the exporter reads *only* the new-era surface -- both the objective
(``objective`` / ``noise_model``, never the retired ``objfunc``, no implicit default;
ADR-0031) **and** the data linkage (``model:`` / ``experiment:`` / ``data:`` /
``condition:`` / ``observable:``; ADR-0028). An ``experiment:`` *is* a PEtab Experiment
(experimentId = the experiment name) carrying its ``data:`` replicates as measurement
rows; a ``condition:`` *is* a PEtab Condition; an ``observable:`` renames a data column
before classification. ``export_job`` **refuses** a legacy (edition 1) job
(:func:`_require_modern_edition`) and a legacy data linkage -- ``model = X : Y.exp`` /
``mutant`` / ``param_scan`` (:func:`_require_new_era_data`) -- rather than reverse-mapping
it. The gate is on the exporter alone; the fitter still runs legacy confs unchanged.
**Conditions/experiments (ADR-0027/0028).** A ``condition:`` referenced by an experiment
becomes a PEtab Condition/Experiment via the surrogate-base ``<p>__REF`` rename of a
fit-and-perturbed parameter (see :mod:`pybnf.petab.conditions`,
:func:`~pybnf.petab.conditions.build_experiment_conditions`); a shared condition emits its
rows once. A dose-response (parameter scan) experiment takes the dual shape (ADR-0046):
each dose becomes a Condition setting the swept parameter + an Experiment, measured at the
scan time (``inf`` for the steady-state default => PEtab time=inf, or a finite ``t_end:``).
The objective and prior surfaces map to PEtab as far as PEtab v2 can express them
(ADR-0023/0031 reversed): the Gaussian/Laplace likelihoods with a ``_SD``-column,
fixed, or column-mean sigma (``chi_sq``/``sos``/``sod``/``ave_norm_sos``), natural-log
Gaussian (``lnnormal`` -> PEtab ``log-normal``), and the
``uniform``/``log-uniform``/``normal``/``laplace`` prior families with their log forms.
A BNGL (``.bngl``) or SBML (``.xml``) model is exported in its own native language
(ADR-0040): a BNGL model PEtab-cleaned, an SBML model carried byte-verbatim with its
observables emitted as ``observableFormula`` expressions from the conf measurement-model
layer (the mirror of the ADR-0036 import). A job may declare **more than one model**
(ADR-0041): each ``experiment:`` names the model it simulates, that model's id is stamped
on the experiment's measurement rows (the ``modelId`` link), free parameters bind across
the union of every model's ids, and ``problem.yaml`` lists every model in its own language
(BNGL + SBML may mix). Everything else raises ``NotImplementedError`` (the boundary is in
code, not silent): an objective PEtab cannot represent (``neg_bin*`` -- removed from v2;
``lognormal`` -- log10 vs PEtab natural log; a free-parameter or relative sigma;
``direct_pass``/``kl``/``wasserstein``); a mean-centred ``lnnormal``, whether the line says
``location = mean`` or the job says ``noise_location = mean`` (#898); the no-prior
``var``/``logvar``; a ``u``-flagged
Uniform, whose box seeds the draw without constraining the search (#736); a ``time_error``
measurement-time marginalization, whose latent sampling time a PEtab measurement row's single
exact ``time`` cannot carry (#738); a fixed-duration equilibration (``equil_t_end:``) on a
model that reads the simulation time, whose clock a PEtab period cannot restart (#896); a
``postprocess`` script (#899); a BNGL action the fit runs ahead of its experiments that can
change their starting state, such as ``setParameter`` (#900); a ``.con``/``.prop``
Constraint; an Antimony (``.ant``) model. The job's ``generate_network`` cap is written into
the exported BNGL model, as the fitter synthesizes it (#901). The
oracle is petab's full ``default_validation_tasks`` via ``Problem.from_yaml`` + the native
``BnglModel`` loader (ADR-0026), wired into the tests; see ADR-0025/0027/0028/0036/0040.
"""
import logging
import math
import re
from dataclasses import dataclass, field
from pathlib import Path
import numpy as np
from .. import edition
from ..data import Data, observed_mean
from ..objective import _NOISE_FAMILIES, _OBJECTIVE_DESUGAR
from ..parameter_record import free_parameter_from_record
from ..parse import ploop
from ..printing import PybnfError
from ..priors import PRIOR_KEYWORD_MAP
from ..pset import (
BNGLModel,
FreeParameter,
INITIALIZATION_PRIOR,
ModelError,
OutOfBoundsException,
)
from ._bngl import parse_model as parse_bngl_model
from ._sbml import parse_model as parse_sbml_model
from .conditions import (
WILDTYPE_CONDITION_ID,
PetabMappingRow,
build_dose_response_conditions,
build_experiment_conditions,
build_preequilibrated_dose_response_conditions,
build_preequilibration_conditions,
is_species_target,
model_time_reads,
species_target_id,
surrogate_name,
write_condition_table,
write_experiment_table,
write_mapping_table,
)
from ._measurement_params import measurement_params_for_replicate, read_measurement_params
from ._tsv import num
from .formula import bngl_body_to_petab_math
from .measurements import (
dose_response_measurement_rows,
measurement_rows_from_data,
write_measurement_table,
)
from .observables import petab_observable_row, write_observable_table
from .parameters import (
EXPORTABLE_PRIOR_KEYWORDS,
petab_parameter_row,
write_parameter_table,
)
logger = logging.getLogger(__name__)
# A noise-model family token (ADR-0031's noise_model grammar / ``_OBJECTIVE_DESUGAR``)
# -> the PEtab v2 noiseDistribution it maps to (ADR-0023 reversed). The LINEAR Gaussian
# / Laplace likelihoods map (``normal`` is the Gaussian alias); the other families are
# explicit boundaries handled in ``_resolve_noise``: ``neg_bin`` was removed from PEtab
# v2, and PyBNF's ``lognormal`` is log10 whereas PEtab's ``log-normal`` is natural log.
# The explicit ``lnnormal`` member is the exact reverse mapping (ADR-0084).
_FAMILY_TOKEN_TO_PETAB_DISTRIBUTION = {
'gaussian': 'normal', 'normal': 'normal', 'lnnormal': 'log-normal',
'laplace': 'laplace'}
# The exported families whose mean and median coincide, so a ``mean`` location scores exactly
# what PEtab's median does (#898). A location-scale family is symmetric on its additive scale,
# so its moment correction vanishes on the LINEAR scale (``ln(base) == 0``): Gaussian's
# ``ln(base)*sigma**2/2`` and Laplace's ``-ln(1 - b**2 t**2)/t`` are both exactly 0 there. On a
# log scale (``lnnormal``) the mean sits above the median and the likelihood changes. Derived
# from the families themselves -- the predicate their own ``supports_profiled_scale`` reads --
# rather than listed, so a family added to the map above lands in the right arm by construction.
_MEAN_IS_MEDIAN_FAMILIES = frozenset(
token for token in _FAMILY_TOKEN_TO_PETAB_DISTRIBUTION
if _NOISE_FAMILIES[token]().additive_on.ln_base == 0.0)
# Free-parameter declaration keywords (the ``(keyword, name)`` tuple keys ``ploop``
# emits). Only ``uniform_var`` exports in chunk 1; the rest raise. The new-era
# ``parameter:`` record is the other declaration spelling and keys on 'parameter' instead
# (ADR-0043); it resolves to one of these same keywords once built (#733).
_VAR_DECL = re.compile(r'(_var$|^var$|^logvar$)')
# The keywords carrying no prior at all (``var``/``logvar``/``lnvar``) -- a flat improper
# prior, which is not a PEtab probability family. Derived from the registry rather than
# listed, so a family added there cannot quietly land in the wrong arm of the refusal
# message; ``has_prior`` is a class attribute, so the family class answers it directly.
_NO_PRIOR_KEYWORDS = frozenset(
keyword for keyword, (family, _scale) in PRIOR_KEYWORD_MAP.items()
if not family.has_prior)
# A legacy ``<name>__FREE`` bind-by-id marker. New-era BNGL binds free parameters by id
# (ADR-0034), so a model carrying this token was not modernized; the exporter refuses it
# rather than ship a ``v1__FREE`` symbol into PEtab (where it would dangle).
_FREE_TOKEN = re.compile(r'\w+__FREE')
# A per-measurement placeholder in a noiseFormula (``noiseParameter1`` / ``observableParameter2``):
# its presence marks a row-varying ``PerMeasurementFormulaSigma`` sigma whose per-row value comes
# from the binding-table sidecar, distinct from a constant ``FormulaSigma`` (ADR-0045). Mirrors
# ``objective._PLACEHOLDER_IN_FORMULA`` / ``import_._PLACEHOLDER``.
_PLACEHOLDER = re.compile(r'(?:observable|noise)Parameter\d')
# ---------------------------------------------------------------------------
# The exporter driver
# ---------------------------------------------------------------------------
[docs]
def export_job(conf_path, out_dir, inline_functions=False):
"""Export the PyBNF job at ``conf_path`` to a PEtab v2 problem in ``out_dir``.
Reads the job's data/conditions/observables from the **new-era surface** (ADR-0028):
a ``model:`` declaration, named ``experiment:`` lines carrying their ``data:`` files
(a PEtab Experiment; multiple files = replicates), ``condition:`` perturbations (a
PEtab Condition), and ``observable:`` column-header overrides. Writes
``parameters.tsv``, ``observables.tsv``, ``measurements.tsv``, a PEtab-clean copy of
the BNGL model, ``problem.yaml``, and -- when some experiment applies a condition --
``conditions.tsv`` / ``experiments.tsv``. Returns the ``out_dir`` path.
``inline_functions`` (default ``False``) is the opt-in expression mode (ADR-0035):
when set, a fitted **function** column emits its body as an ``observableFormula``
expression (translated to PEtab math, requires the ``pybnf[petab]`` extra) instead of
the bare model name, producing a model-portable problem and the round-trip oracle the
importer's synthesis is graded against. The default stays bare-name, lossless,
byte-stable, and ``petab``-free; an observable column is never inlined.
New-era only (ADR-0028 Chunk 5c, "refuse legacy everything"): a job that binds data
the legacy way (``model = X : Y.exp`` / ``mutant`` / ``param_scan``) is refused, as is
a legacy-edition job (:func:`_require_modern_edition`). Raises ``NotImplementedError``
at every documented boundary and ``PybnfError`` for a malformed/unsupported job (see
the module docstring).
"""
conf_path = Path(conf_path)
out_dir = Path(out_dir)
out_dir.mkdir(parents=True, exist_ok=True)
conf = _read_conf_dict(conf_path)
_require_modern_edition(conf)
models = _resolve_models(conf)
languages = {mf: _model_language(mf) for mf in models}
_require_new_era_data(conf, models)
noise = _resolve_noise(conf)
per_obs_noise = _resolve_per_observable_noise(conf)
_reject_cumulative(conf)
_reject_time_error(conf)
_reject_normalization(conf)
_reject_postprocess(conf)
free_params = _free_parameters_from_conf(conf)
# An estimated (`fit`) sigma exports as a bare-id noiseFormula naming an estimated PEtab
# parameter (#439), so its noise scale must be a DECLARED free parameter (with bounds/prior
# to write). The legacy whole-fit `chi_sq_dynamic` synthesizes an IMPLICIT `sigma__FREE` the
# user never declared, so it has nothing to write -- still a boundary.
_declared_free = {fp.name for fp in free_params}
for _dist, _verb, _arg in [noise, *per_obs_noise.values()]:
if _verb == 'fit' and _arg not in _declared_free:
raise NotImplementedError(
f"The objective's estimated sigma '{_arg}' is an implicit noise parameter "
f"(e.g. chi_sq_dynamic's sigma__FREE) that is not a declared free parameter, so "
f"it has no bounds/prior to export as a PEtab estimated parameter. Declare it as "
f"a free parameter (as a per-observable `fit` noise scale does) to export it. "
f"ADR-0021/0023, #439.")
# A registry of per-language model views (ADR-0040/0041): a job may mix BNGL + SBML,
# each read once and threaded through the language-agnostic classification below.
registry = {mf: _read_model(mf, conf_path.parent / mf, languages[mf]) for mf in models}
# Observation-layer nuisance free parameters (ADR-0034/0044/0045): a measurement scale, a
# noise coefficient, or a row-varying per-row sigma -- a free parameter referenced by a
# measurement-model / noise formula or a binding-table token but NOT a model entity. The
# bind-by-id check admits them as unbound estimated parameters (else they read as typos);
# the rest must still bind to a model id.
model_ids = set().union(*(set(v.parameters) for v in registry.values()))
nuisances = (_referenced_nuisance_symbols(conf, conf_path, noise, per_obs_noise)
& {fp.name for fp in free_params}) - model_ids
free_to_model = _resolve_free_to_model(free_params, registry, models, nuisances)
fit_model_params = set(free_to_model.values()) & model_ids
(observable_rows, measurement_rows, condition_rows, experiment_rows,
surrogate_params, mapping_rows) = _export_new_era(
conf, conf_path, models, registry, noise, per_obs_noise, fit_model_params,
inline_functions)
parameter_rows = _parameter_rows(
free_params, free_to_model, surrogate_params, registry, models)
# Each model is emitted in its own native language (ADR-0040): a BNGL model is
# PEtab-cleaned (its actions reduced to the network definition the fit used, #485/#900/
# #901); an SBML model is carried byte-verbatim (the measurement model lives in the
# observables table, never a model-file edit -- ADR-0036). Cleaned before any file is
# written, so a refused actions block leaves no half-written problem behind.
methods = _experiment_methods_by_model(conf, models)
model_texts = {
mf: (clean_model_for_petab(
registry[mf].text, Path(mf).name,
generate_network_options=conf.get('generate_network'),
experiment_methods=methods[mf])
if languages[mf] == 'bngl' else registry[mf].text)
for mf in models}
write_parameter_table(parameter_rows, out_dir / 'parameters.tsv')
write_observable_table(observable_rows, out_dir / 'observables.tsv')
write_measurement_table(measurement_rows, out_dir / 'measurements.tsv')
if condition_rows:
write_condition_table(condition_rows, out_dir / 'conditions.tsv')
if experiment_rows:
write_experiment_table(experiment_rows, out_dir / 'experiments.tsv')
if mapping_rows:
write_mapping_table(mapping_rows, out_dir / 'mapping.tsv')
for mf in models:
(out_dir / Path(mf).name).write_text(model_texts[mf])
# One model_files entry per model (ADR-0041), modelId = the file stem, in declaration
# order (so a re-export reproduces the same problem.yaml model_files ordering).
model_yaml = [(Path(mf).stem, Path(mf).name, languages[mf]) for mf in models]
write_problem_yaml(out_dir / 'problem.yaml', model_yaml,
has_conditions=bool(condition_rows),
has_experiments=bool(experiment_rows),
has_mapping=bool(mapping_rows))
return out_dir
# ---------------------------------------------------------------------------
# New-era surface reading (ADR-0028): export becomes transcription
# ---------------------------------------------------------------------------
def _has_new_era_data(conf):
"""True iff the job binds data via the new-era ``experiment:`` surface (ADR-0028).
A fully new-era conf introduces data through ``('experiment', name)`` entries (a PEtab
Experiment carrying its ``data:``), never the legacy ``model = X : Y.exp`` linkage.
"""
return any(isinstance(k, tuple) and len(k) == 2 and k[0] == 'experiment'
for k in conf)
def _require_new_era_data(conf, models):
"""Refuse a legacy data linkage -- the PEtab v2 exporter reads only the new-era data
surface (ADR-0028 Chunk 5c, "refuse legacy everything").
The exporter is already new-era-gated on the *objective* (``_require_modern_edition``);
now that ADR-0028's data surface exists, the data linkage is held to the same standard.
A job must introduce data through named ``experiment:`` lines (a PEtab Experiment
carrying its ``data:``); the legacy filename->suffix binding (``model = X : Y.exp``),
``mutant`` lines, and ``param_scan`` actions are refused rather than silently read.
Mixing the two (a ``mutant``/``param_scan`` line, or data on any ``model =`` line,
alongside ``experiment:``) is refused too, so a legacy line is never silently dropped.
The gate is on the *exporter* only; the fitter still runs legacy confs unchanged.
"""
legacy_data = any(conf.get(mf) for mf in models) # model = X : Y.exp data list
legacy_features = [k for k in ('mutant', 'param_scan') if k in conf]
if not _has_new_era_data(conf):
raise NotImplementedError(
"The PEtab v2 exporter reads the new-era data surface; this job binds data "
"the legacy way (model = X : Y.exp / mutant / param_scan), which is refused "
"(ADR-0028, 'refuse legacy everything'). Re-author it on the new-era surface: "
"declare the model with 'model:', bind data through a named 'experiment:' (one "
"'data:' file per replicate), and write perturbations as 'condition:' lines. "
"The fitter still runs the legacy form; only export requires the new surface.")
if legacy_data or legacy_features:
legacy = (['data on the legacy model = X : Y.exp line'] if legacy_data else []) \
+ [f"a '{k}' line" for k in legacy_features]
raise NotImplementedError(
f"This job mixes the new-era 'experiment:' surface with legacy data linkage "
f"({', '.join(legacy)}). Use the new-era surface exclusively -- move all data "
f"into 'experiment:'/'data:' and all perturbations into 'condition:' -- so no "
f"legacy line is silently ignored on export (ADR-0028, #423).")
def _export_new_era(conf, conf_path, models, registry, noise, per_obs_noise,
fit_model_params, inline_functions=False):
"""Read a job's data/conditions/observables from the **new-era surface** (ADR-0028).
Export is *transcription*: an ``experiment:`` is a PEtab Experiment (experimentId =
the experiment name) carrying its ``data:`` replicates as measurement rows; an
``observable:`` line renames a data column to a model entity before classification.
``models`` is the ordered list of model files and ``registry`` their per-language views
(ADR-0041): each experiment names the model it simulates, that model's id is stamped on
its measurement rows' ``modelId`` (omitted when the job is single-model), and a column
is classified against its experiment's model. ``noise`` is the whole-fit base and
``per_obs_noise`` the ``{column: (dist, verb, arg)}`` per-observable overrides
(ADR-0021/0045): a column with an override takes its own sigma source, the rest the base.
``inline_functions`` is threaded to :func:`_observable_rows` (ADR-0035 inlining). Returns
``(observable_rows, measurement_rows, condition_rows, experiment_rows, surrogate_params)``.
A ``condition:`` referenced by an experiment becomes a PEtab Condition (the
surrogate-base machinery of ADR-0027, generalized by
:func:`~pybnf.petab.conditions.build_experiment_conditions`): a fit-and-perturbed
parameter is renamed to ``<p>__REF`` in the parameter table and pinned in every
experiment's Condition. A **parameter-scan** (dose-response) experiment takes the
dual shape (ADR-0046): each dose of its dose axis (the sorted union of its replicates'
doses, #895) becomes a Condition setting the swept parameter (and re-pinning M, #892) + an
Experiment, measured at the scan time (``inf`` for the
steady-state default => PEtab time=inf, or a finite ``t_end:``), via
:func:`~pybnf.petab.conditions.build_dose_response_conditions` +
:func:`~pybnf.petab.measurements.dose_response_measurement_rows`. With no referenced
conditions the surrogate set is empty, so a single wildtype time course is byte-identical
to the chunk-1 base.
A **pre-equilibration** experiment (``preequilibrate:``, ADR-0052) takes a third shape: a
two-period Experiment (a ``time = -inf`` steady-state period under the pre-equilibration
condition -- or a ``time = -T`` period for a fixed-duration ``equil_t_end: T``, #896 -- + a
``time = 0`` period under the measurement condition), built by
:func:`~pybnf.petab.conditions.build_preequilibration_conditions`. Its measurements are
tagged exactly like a time course's (the data grid at times >= 0; the equilibration period
carries no measurements).
A **pre-equilibrated dose-response** (``preequilibrate:`` + ``parameter_scan``, ADR-0062 --
the preincubate -> wash -> dose-scan protocol) takes a fourth shape: N two-period Experiments,
one per dose, whose ``time = 0`` measurement period applies BOTH a shared wash condition and a
per-dose swept-parameter condition, built by
:func:`~pybnf.petab.conditions.build_preequilibrated_dose_response_conditions`. Its
measurements pivot exactly like a plain dose-response's (tagged ``<stem>_<i>`` at the scan
time). A species ``setConcentration`` target in any condition (a wash/bolus, ADR-0062) is
aliased through the **mapping table**: a BNGL species pattern is not a valid PEtab id, so
:func:`_species_id_map` synthesizes a ``species_<...>`` id (``petabEntityId``) for it, the
condition targets that id, and ``mapping_rows`` (``petab_id -> pattern``) is emitted as
``mapping.tsv``. Returns ``mapping_rows`` as its sixth element (empty for a job with no
species condition, so ``mapping.tsv`` stays absent).
"""
experiments = _read_experiments(conf, conf_path, models)
overrides = _read_observable_overrides(conf)
all_datas = [d for exp in experiments for d in exp['datas']]
_apply_observable_overrides(all_datas, overrides)
measurement_models = _read_measurement_models(conf)
observable_rows, column_to_observable_id, column_means = _observable_rows(
experiments, registry, noise, per_obs_noise, inline_functions, measurement_models)
def _noise_values_for(exp):
# A column-mean sigma that differs between experiments (#894) is written row by row:
# every measurement row of this experiment carries this experiment's own mean.
return {col: means[exp['name']] for col, means in column_means.items()
if exp['name'] in means}
# Four PEtab experiment shapes (ADR-0046/0052/0062): a time course is one Experiment over a
# referenced Condition; a dose-response (parameter_scan) is N Conditions (each sets the swept
# parameter to one dose) + N Experiments measured at the scan time (inf for steady state); a
# pre-equilibration experiment is a two-period Experiment (a -inf steady-state period + a
# time=0 measurement period -- ADR-0052); a PRE-EQUILIBRATED dose-response (a preincubate ->
# wash -> dose-scan protocol, ADR-0062) is N two-period Experiments, one per dose, whose
# measurement period applies both a shared wash condition and a per-dose condition. They build
# independently and concatenate; a `preequilibrate:` field splits an experiment off the plain
# time-course/dose-response buckets, and its `type` splits pre-equilibration (time course) from
# a pre-equilibrated scan.
pe_experiments = [exp for exp in experiments
if exp['preequilibrate'] is not None and exp['type'] == 'time_course']
pdr_experiments = [exp for exp in experiments
if exp['preequilibrate'] is not None and exp['type'] == 'parameter_scan']
tc_experiments = [exp for exp in experiments
if exp['type'] == 'time_course' and exp['preequilibrate'] is None]
dr_experiments = [exp for exp in experiments
if exp['type'] == 'parameter_scan' and exp['preequilibrate'] is None]
_refuse_fixed_equilibration_of_time_dependent_models(pe_experiments + pdr_experiments,
registry)
conditions, condition_models = _read_conditions(conf, models, registry)
referenced = {exp['condition'] for exp in tc_experiments if exp['condition'] is not None}
# A pre-equilibration / pre-equilibrated-scan experiment references its pre-equilibration
# condition AND (optionally) its measurement (wash) condition NOT via the time-course
# ``condition:`` path, so add both to ``referenced`` -- else _read_conditions drops them as
# "unused" (ADR-0052/0062).
for exp in pe_experiments + pdr_experiments:
referenced.add(exp['preequilibrate'])
if exp['condition'] is not None:
referenced.add(exp['condition'])
undefined = referenced - set(conditions)
if undefined:
raise PybnfError(
f"Experiment(s) reference undefined condition(s) {sorted(undefined)}; define "
f"each with a 'condition:' line.")
# A condition named 'wildtype' would be written as conditionId cond_wildtype, the id the
# exporter reserves for its synthesized base condition (ADR-0027). A reader cannot tell the
# two apart by id -- PyBNF 1.8.1's importer dropped every cond_wildtype row (#905) -- and the
# builders write each id once, so where the base is also needed (M non-empty: a wildtype time
# course, a wash-out, a `none` equilibration period, #906) one set of rows would silently serve
# both. So the name is refused whenever an experiment applies such a condition, not only when
# the base is emitted -- a `perturbations: none` one too, measured or pre-equilibration, so the
# rule is simply that the name is reserved.
applied = referenced | {exp['unperturbed_condition'] for exp in experiments
if exp['unperturbed_condition'] is not None}
clash = sorted(c for c in applied if f'cond_{c}' == WILDTYPE_CONDITION_ID)
if clash:
raise PybnfError(
f"Condition '{clash[0]}' cannot be exported to PEtab: it would be written as "
f"conditionId '{WILDTYPE_CONDITION_ID}', which the exporter reserves for the base "
f"condition it synthesizes for wildtype and wash-out experiments, and a PEtab reader "
f"may take its targets for that base and drop them. Rename the condition (any other "
f"name is written as 'cond_<name>').")
# A condition belongs to one model (ADR-0041 addendum), and the fitter looks an experiment's
# conditions up on the experiment's OWN model only (config.py::_resolve_experiment_data_key /
# _preequilibration_perturbations), so an experiment applying another model's condition is a
# job the fitter refuses. Refuse it here too: the condition's model is what its fixed-target
# relative ops are computed against (#897), so there is no single right base otherwise.
for exp in experiments:
for c in (exp['condition'], exp['preequilibrate'], exp['unperturbed_condition']):
if c is not None and c in condition_models and condition_models[c] != exp['model']:
raise PybnfError(
f"Experiment '{exp['name']}' simulates model '{exp['model']}' but applies "
f"condition '{c}', which belongs to model '{condition_models[c]}'. A "
f"condition perturbs only the model it names (ADR-0041); declare the "
f"condition for '{exp['model']}' or apply it to an experiment on that model.")
nominal_of = _condition_nominal_of(registry, condition_models)
# A species-target condition (setConcentration -- a wash/bolus, ADR-0062) exports to a PEtab v2
# condition whose target is a species *amount*: a BNGL pattern is not a valid PEtab id, so each
# referenced species pattern is aliased to a synthesized ``species_<...>`` id via the mapping
# table (petabEntityId -> the pattern), and the condition targets that id. Built once, globally.
species_id_of = _species_id_map(
var for name in referenced for var, _op, _val in conditions.get(name, [])
if is_species_target(var))
# A species setConcentration is inline-only within a pre-equilibration protocol (ADR-0062, the
# fitter's rule): it perturbs a species amount mid-protocol, not a mutant parameter block. A
# plain time-course experiment applying a species condition has no such inline phase, so refuse
# it with a clear message (the fitter rejects it too) rather than mis-emit.
tc_species_conds = sorted(
exp['condition'] for exp in tc_experiments
if exp['condition'] is not None
and any(is_species_target(var) for var, _o, _v in conditions.get(exp['condition'], [])))
if tc_species_conds:
raise NotImplementedError(
f"Condition(s) {tc_species_conds} set a species amount (setConcentration) but are "
f"applied to a plain time-course experiment. A species setConcentration is inline-only "
f"within a pre-equilibration protocol (ADR-0062); apply it through a 'preequilibrate:' "
f"experiment (as an incubate or a wash), not a bare 'condition:' time course.")
unused = set(conditions) - referenced
if unused:
# An unused condition emits no PEtab rows (the fitter would not apply it either);
# skip it with a debug log rather than warning (ADR-0028 Chunk 5, decision 4).
logger.debug("Conditions defined but referenced by no experiment (skipped): %s",
sorted(unused))
# The surrogate set M (ADR-0027) is problem-global: removing a fit-and-perturbed param from
# the parameter table makes its model name a pure condition target, so EVERY experiment of
# EVERY shape must set it before its simulation starts -- the time-course, wildtype,
# pre-equilibration, dose-response and pre-equilibrated dose-response conditions alike (#443,
# #892). A pre-equilibration or pre-equilibrated-scan condition that perturbs a fit param
# therefore contributes to M too; thread that contribution into the time-course builder's M
# via extra_surrogate, so every builder below pins the same M.
pe_surrogate = {
var for exp in pe_experiments + pdr_experiments
for c in ([exp['preequilibrate']]
+ ([exp['condition']] if exp['condition'] is not None else []))
for var, _op, _val in conditions[c] if var in fit_model_params}
condition_rows, experiment_rows, surrogate_params, experiment_to_id = \
build_experiment_conditions(
[(exp['name'], exp['condition']) for exp in tc_experiments],
conditions, fit_model_params, nominal_of,
extra_surrogate=pe_surrogate)
# Pre-equilibration experiments -> two-period Experiments (ADR-0052): a -inf steady-state
# period (or a -T period for a fixed equil_t_end: T, #896) under the pre-equilibration
# condition + a time=0 period under the measurement condition. They share the
# problem-global M (every period re-pins M -- #443): a fit-parameter perturbation in a
# pre-equilibration period emits its surrogate op and every other period re-pins the base;
# a wash-out re-pins M via the synthesized cond_wildtype base.
# existing_condition_ids dedups a condition shared with a time course and the wildtype base.
if pe_experiments:
pe_condition_rows, pe_experiment_rows, pe_experiment_to_id = \
build_preequilibration_conditions(
[(exp['name'], exp['preequilibrate'], exp['condition'], exp['equil_t_end'])
for exp in pe_experiments],
conditions, nominal_of,
surrogate=surrogate_params,
existing_condition_ids={r.condition_id for r in condition_rows},
species_id_of=species_id_of)
condition_rows += pe_condition_rows
experiment_rows += pe_experiment_rows
experiment_to_id.update(pe_experiment_to_id)
# Pre-equilibrated dose-response experiments (ADR-0062): N two-period Experiments per scan, a
# -inf (or -equil_t_end, #896) pre-equilibration period + a measurement period applying the
# shared wash condition and a per-dose swept-parameter condition. They share the
# problem-global M (#892): the pre-equilibration condition sets all of M in its leading
# period, and PEtab v2, like the fitter, carries those values into the measurement period. A
# wash condition re-pins M as every other condition does, and _refuse_wash_re_pins refuses
# the cases where that re-pin is wrong.
if pdr_experiments:
_refuse_wash_re_pins(pdr_experiments, conditions, surrogate_params)
pdr_condition_rows, pdr_experiment_rows, pdr_ids_by_name = \
build_preequilibrated_dose_response_conditions(
[(exp['name'], exp['preequilibrate'], exp['condition'],
_swept_param(exp), _dose_axis(exp), exp['scan_time'], exp['equil_t_end'])
for exp in pdr_experiments],
conditions, nominal_of,
species_id_of=species_id_of,
existing_condition_ids={r.condition_id for r in condition_rows},
surrogate=surrogate_params)
condition_rows += pdr_condition_rows
experiment_rows += pdr_experiment_rows
# Per-point numeric noiseParameters are emitted only when a column's sigma comes from a
# data column (the read_exp_file placeholder source) or is a column mean that differs
# between experiments (_noise_values_for above, #894); a fixed / single-mean column-mean /
# formula sigma is carried inline in noiseFormula, so the measurement export must not read
# _SD then (it would leave a noiseParameters override with no placeholder to bind to).
# With per-observable overrides the suffix is **per column** (ADR-0045): each column uses
# its own sigma source (its override, else the whole-fit base) to decide whether it reads
# a _SD companion.
def _sd_suffix_for(col):
_dist, verb, arg = per_obs_noise.get(col, noise)
return arg if verb == 'read_exp_file' else None
sd_suffix = {col: _sd_suffix_for(col) for col in column_to_observable_id}
# The modelId link (ADR-0041): each experiment stamps its model's stem onto its
# measurement rows. Single-model -> '' (the column is dropped on write, byte-stable).
multi_model = len(models) > 1
measurement_rows = []
# A pre-equilibration experiment's measurements are tagged exactly like a time course's
# (the data grid at times >= 0 under its experimentId): the -inf equilibration period
# carries no measurements, and PEtab resolves the data times into the time=0 period (ADR-0052).
for exp in tc_experiments + pe_experiments:
eid = experiment_to_id[exp['name']]
model_id = Path(exp['model']).stem if multi_model else ''
# Each replicate Data contributes its own rows under the one experiment (PEtab
# models replicates as repeated rows -- no need to pre-stack as config.py does). A
# row-varying placeholder's per-row token comes from that replicate's slice of the
# experiment's measurement_params sidecar (ADR-0083); a legacy four-column sidecar is
# selected for every replicate and therefore keeps its original shared-token semantics.
for replicate, data in enumerate(exp['datas']):
cmap = {c: o for c, o in column_to_observable_id.items() if c in data.cols}
measurement_rows += measurement_rows_from_data(
data, cmap, experiment_id=eid, sd_suffix=sd_suffix, model_id=model_id,
measurement_params=measurement_params_for_replicate(
exp['measurement_params'], replicate),
noise_values=_noise_values_for(exp))
# Dose-response (ADR-0046): each dose of the experiment's dose axis -- the sorted union of its
# replicates' doses, the grid the fitter scans (#895) -- becomes its own Condition (setting the
# swept parameter, and re-pinning M, #892) + Experiment, and the observable columns become
# measurements at the scan time (inf => steady state). Each data row is tagged with the
# experiment of its OWN dose, never by its row position, so a replicate whose doses are
# reordered, missing or extra pairs every measurement with the dose it was taken at, as the
# fitter does. The swept-parameter column is the scan axis, not a measurement, so it is
# dropped from the column map.
for exp in dr_experiments:
stem = exp['name']
model_id = Path(exp['model']).stem if multi_model else ''
scan_time = exp['scan_time']
swept_param = _swept_param(exp)
dose_axis = _dose_axis(exp)
dr_conditions, dr_experiment_rows, experiment_ids = build_dose_response_conditions(
stem, swept_param, dose_axis, scan_time, surrogate=surrogate_params)
condition_rows += dr_conditions
experiment_rows += dr_experiment_rows
experiment_id_of_dose = dict(zip(dose_axis, experiment_ids))
for data in exp['datas']:
cmap = {c: o for c, o in column_to_observable_id.items()
if c in data.cols and c != swept_param}
# A column-mean sigma is the whole scan's mean, not a per-dose one: the scan is
# ONE PyBNF experiment however many PEtab experiments its doses become (#894).
measurement_rows += dose_response_measurement_rows(
data, cmap, swept_param, experiment_id_of_dose, scan_time,
sd_suffix=sd_suffix, model_id=model_id,
noise_values=_noise_values_for(exp))
# Pre-equilibrated dose-response measurements (ADR-0062): tagged <stem>_<i> at the scan time,
# exactly like a plain dose-response, so the same pivot applies -- each row by its own dose
# (#895). The per-experiment experiment ids come from the builder above, aligned with the same
# dose axis; the swept-parameter column is the scan axis, not a measurement.
for exp in pdr_experiments:
model_id = Path(exp['model']).stem if multi_model else ''
scan_time = exp['scan_time']
swept_param = _swept_param(exp)
experiment_id_of_dose = dict(zip(_dose_axis(exp), pdr_ids_by_name[exp['name']]))
for data in exp['datas']:
cmap = {c: o for c, o in column_to_observable_id.items()
if c in data.cols and c != swept_param}
measurement_rows += dose_response_measurement_rows(
data, cmap, swept_param, experiment_id_of_dose, scan_time,
sd_suffix=sd_suffix, model_id=model_id,
noise_values=_noise_values_for(exp))
# The species-amount mapping table (ADR-0062): one row per referenced species pattern, in
# first-appearance order (petabEntityId -> the BNGL pattern). Empty for a job with no species
# setConcentration condition, so mapping.tsv / problem.yaml mapping_files stay absent then.
mapping_rows = [PetabMappingRow(pid, pattern) for pattern, pid in species_id_of.items()]
return observable_rows, measurement_rows, condition_rows, experiment_rows, \
surrogate_params, mapping_rows
def _read_experiments(conf, conf_path, models):
"""Read + resolve the new-era ``experiment:`` entries from the raw ``ploop`` dict.
Each ``('experiment', name)`` entry is ``{'data': [files], 'condition': c?, 'model':
mf?, 'type': t?, 'method': m?, 't_end': t?, 'preequilibrate': p?, 'measurement_params':
mp?, 'equil_t_end': T?}``. ``models`` is the ordered list of the job's model files. Returns a
list (declaration order) of dicts ``{'name', 'condition', 'unperturbed_condition',
'model': model_file, 'datas': [Data, ...], 'data_files': [str, ...], 'type', 'scan_time',
'preequilibrate': cond?, 'measurement_params': table?, 'equil_t_end': T?}`` (``T`` the fixed
equilibration duration, :func:`_equil_t_end`; a measured ``perturbations: none`` condition
reads as ``'condition': None`` with its name in ``'unperturbed_condition'``, #906) -- the
``data:`` files (``data_files``, as written in the conf, for
error messages) read as individual
:class:`~pybnf.data.Data` replicates (PEtab models replicates as repeated measurement
rows, so they are not pre-stacked), each experiment's resolved model
(:func:`_resolve_experiment_model`, ADR-0041), the inferred ``type`` (``'time_course'`` or
``'parameter_scan'``), the dose-response ``scan_time`` (``inf`` for the steady-state
default, a finite ``t_end:`` otherwise; ``None`` for a time course -- ADR-0046), and its
row-varying per-measurement binding table read from the ``measurement_params:`` sidecar
(``{column: {placeholder: {key: token}}}``, where ``key`` is ``(replicate, time)`` in the
ADR-0083 format or a shared bare time in the legacy format; ``None`` when absent). Raises:
* the ambiguous-model error if an experiment names no model but the job has more than
one (mirrors ``config.py::_resolve_experiment_model``);
* a not-yet-supported boundary for a parameter_scan that also names a ``condition:``
(a dose-response already makes each dose its own condition -- ADR-0046), or a
``preequilibrate:`` parameter_scan (a scan after equilibration has no export route);
* a constraint refusal for non-``.exp`` data: BPSL ``.con``/``.prop`` constraints are
PyBNF-native with no core-PEtab representation, so an experiment carrying them cannot
be exported (the fitter still runs it -- ADR-0028 addendum).
"""
stem_to_model = {Path(mf).stem: mf for mf in models}
# A `perturbations: none` condition (#906, ADR-0150) changes nothing, so as the measured
# `condition:` it is exactly an omitted one and exports as such (read here as None). As
# `preequilibrate:` it keeps its name -- that is what makes the experiment two-period -- and
# the builders write its -inf period as the model as is.
unperturbed = {k[1] for k, v in conf.items()
if isinstance(k, tuple) and len(k) == 2 and k[0] == 'condition' and not v[1]}
experiments = []
for key, fields in conf.items():
if not (isinstance(key, tuple) and len(key) == 2 and key[0] == 'experiment'):
continue
name = key[1]
condition = fields.get('condition')
# Kept, so the model-ownership check can still refuse another model's `none` condition
# exactly as the fitter does (which looks it up on the experiment's own model).
unperturbed_condition = None
if condition in unperturbed:
condition, unperturbed_condition = None, condition
model_file = _resolve_experiment_model(name, fields.get('model'), models,
stem_to_model)
data_files = fields.get('data', [])
if not data_files:
raise PybnfError(f"Experiment '{name}' declares no 'data:' files.")
non_exp = [f for f in data_files if not f.endswith('.exp')]
if non_exp:
raise NotImplementedError(
f"Experiment '{name}' carries BPSL constraint data ({non_exp}). PyBNF "
f"constraints (.con/.prop) are a native qualitative-fitting feature with no "
f"core-PEtab v2 representation, so this experiment cannot be exported "
f"(ADR-0028). The fitter still runs it; only export is refused. Drop the "
f"constraint file(s) from 'data:' to export the quantitative measurements "
f"alone.")
datas = [Data(file_name=str(conf_path.parent / f)) for f in data_files]
exp_type = _experiment_type(name, datas[0], fields.get('type'))
preequilibrate = fields.get('preequilibrate')
# A pre-equilibration experiment (ADR-0052) is measured AFTER an unmeasured equilibration
# phase -> a PEtab two-period Experiment. A time-course measured phase gives a single
# measurement period; a parameter_scan measured phase (a preincubate->wash->dose-scan
# protocol, ADR-0062/#477) gives N two-period experiments, one per dose, with a
# multi-condition measurement period -- both export routes below.
# A parameter_scan (dose-response) experiment's measurement time is its scan endpoint
# (ADR-0046): inf for the steady-state default (PEtab time=inf), or a finite ``t_end:``.
# A time course derives its grid from the data, so ``t_end:`` is inert there.
scan_time = None
if exp_type == 'parameter_scan':
t_end = fields.get('t_end')
scan_time = float(t_end) if t_end is not None else float('inf')
# A PLAIN dose-response makes each dose its own condition (ADR-0046), so a named
# condition on it has no export route. A PRE-EQUILIBRATED dose-response (ADR-0062),
# by contrast, names its measurement (wash) condition, applied alongside the per-dose
# condition in the measurement period -- that route exists, so allow it there.
if condition is not None and preequilibrate is None:
raise NotImplementedError(
f"Experiment '{name}' is a parameter_scan that also names a condition "
f"('{condition}'). A dose-response already makes each dose its "
f"own condition (ADR-0046); combining it with a named condition has no "
f"export route yet.")
measurement_params = None
mp_file = fields.get('measurement_params')
if mp_file:
measurement_params = read_measurement_params(conf_path.parent / mp_file)
experiments.append({'name': name, 'condition': condition,
'unperturbed_condition': unperturbed_condition,
'model': model_file, 'datas': datas, 'data_files': list(data_files),
'type': exp_type,
'scan_time': scan_time, 'preequilibrate': preequilibrate,
'measurement_params': measurement_params,
'equil_t_end': _equil_t_end(name, fields, preequilibrate)})
return experiments
def _experiment_methods_by_model(conf, models):
"""``{model_file: {experiment_name: method}}`` -- each experiment's simulation ``method:``
(``ode`` when unset, as ``config.py`` defaults it), grouped by the model it simulates.
PEtab has no simulation method, so the tables never read this; the model cleaner does.
The fitter builds each experiment's simulation differently by method: a network-free
(``nf``) one is emitted with no ``resetConcentrations()`` before it and does not make the
fitter generate a network (``BNGLModel.add_action``). So whether a hand-written
simulation action can change what an experiment starts from (#900), and whether the fit
generates a network at all (#901), depend on it."""
stem_to_model = {Path(mf).stem: mf for mf in models}
methods = {mf: {} for mf in models}
for key, fields in conf.items():
if isinstance(key, tuple) and len(key) == 2 and key[0] == 'experiment':
mf = _resolve_experiment_model(key[1], fields.get('model'), models, stem_to_model)
methods[mf][key[1]] = str(fields.get('method', 'ode')).lower()
return methods
def _equil_t_end(name, fields, preequilibrate):
"""An experiment's fixed equilibration duration (``equil_t_end:``) as a float, or ``None``.
The fitter runs a ``preequilibrate:`` experiment's unmeasured phase for exactly this long
instead of to steady state (``config.py::_build_preequilibration_action`` ->
``pset.py::_append_preequilibration_actions``), so it exports as a leading PEtab period at
``time = -equil_t_end`` rather than ``-inf`` (#896). That period is distinct from the
measured one at 0 only for a finite, positive duration: a zero duration would start the two
periods at the same time (PEtab would apply both conditions at once), and a negative or
infinite one is no duration at all -- each is refused with the experiment named. Without
``preequilibrate:`` the fitter never reads the field, so neither does the export."""
raw = fields.get('equil_t_end')
if raw is None or preequilibrate is None:
return None
t = float(raw)
if not (math.isfinite(t) and t > 0):
raise PybnfError(
f"Experiment '{name}' sets equil_t_end: {raw}, but a fixed equilibration duration must "
f"be a finite positive time. Omit equil_t_end to equilibrate to steady state, or give "
f"the duration of the unmeasured equilibration phase.")
return t
def _refuse_fixed_equilibration_of_time_dependent_models(experiments, registry):
"""Refuse a fixed-duration equilibration (``equil_t_end: T``) on a model that reads time.
The exported leading period runs on ``t`` in ``[-T, 0]`` (PEtab v2: a period lasts from its
start time to the next period's), but the fitter runs the same phase on ``[0, T]`` and
restarts the clock at 0 for the measured phase. For an autonomous model the two are the same
simulation; for one whose rates, functions, or events read the time they are not, and PEtab
v2 has no way to restart the clock between periods, so there is no exact representation
(#896). ``experiments`` are the pre-equilibration / pre-equilibrated-scan experiment dicts."""
for exp in experiments:
if exp['equil_t_end'] is None:
continue
model_file = exp['model']
reads = model_time_reads(registry[model_file].text, _model_language(model_file))
if reads:
raise NotImplementedError(
f"Experiment '{exp['name']}' equilibrates for a fixed duration "
f"(equil_t_end: {num(exp['equil_t_end'])}) on model '{model_file}', which reads "
f"the simulation time ({'; '.join(reads)}). PyBNF runs that equilibration from "
f"t = 0 to t = {num(exp['equil_t_end'])} and restarts the clock at 0 for the "
f"measured phase, but a PEtab v2 equilibration period runs from "
f"t = -{num(exp['equil_t_end'])} to 0, so the model sees different times and the "
f"exported problem would not reproduce the fit. Run the job natively, or export "
f"it with a model that does not read time (#896).")
def _resolve_experiment_model(name, ref, models, stem_to_model):
"""The model file an experiment simulates (ADR-0041), the export peer of
``config.py::_resolve_experiment_model``.
With an explicit ``model:`` ref, resolve it by filename stem (the ``model_files`` key);
an unknown ref is a typo. With no ref, default to the sole model when the job declares
exactly one; under more than one model an unnamed experiment is ambiguous -- the
exporter requires the ``model:`` field rather than guessing which model produced the
data (matching the fitter's rule).
"""
if ref is not None:
stem = Path(ref).stem
if stem not in stem_to_model:
raise PybnfError(
f"Experiment '{name}' names model '{ref}', but the job declares no model "
f"with id '{stem}' (declared model ids: {sorted(stem_to_model)}).")
return stem_to_model[stem]
if len(models) == 1:
return models[0]
raise PybnfError(
f"Experiment '{name}' does not name a model, but the job declares {len(models)} "
f"models ({models}). Add 'model: <file>' to the experiment to say which model it "
f"simulates (ADR-0041).")
def _experiment_type(name, data, explicit_type):
"""Infer ``'time_course'`` vs ``'parameter_scan'`` from a ``Data``'s independent
variable (``time`` => time_course; otherwise the indvar names a swept parameter =>
parameter_scan), unless ``type:`` states it. Mirrors
``config.py::_infer_experiment_type`` (a scan exports each dose as a steady-state
Condition/Experiment -- ADR-0046).
A ``steady_state`` experiment (ADR-0086, #521) exports on the **time-course route**:
its measurement time already IS ``inf`` in the ``.exp``, which is exactly PEtab's
steady-state time, so the rows need no special casing -- only the ``type:`` token has to
be accepted here (the fitter's inference and the exporter's differ only in that the
fitter must pick a simulation for it)."""
if explicit_type is not None:
t = explicit_type.lower()
if t in ('time_course', 'timecourse', 'steady_state', 'steadystate'):
return 'time_course'
if t in ('parameter_scan', 'param_scan', 'parameterscan'):
return 'parameter_scan'
raise PybnfError(
f"Experiment '{name}' has unrecognized type '{explicit_type}' (use "
f"'time_course').")
indvar = data.indvar if data.indvar is not None else _independent_variable(data)
return 'time_course' if indvar.lower() == 'time' else 'parameter_scan'
def _read_observable_overrides(conf):
"""The new-era ``observable: <entity>, column: <header>`` overrides as
``{entity: header}`` (ADR-0028 Chunk 4) -- the renames applied before classification."""
return {k[1]: v for k, v in conf.items()
if isinstance(k, tuple) and len(k) == 2 and k[0] == 'observable'}
def _read_measurement_models(conf):
"""The new-era ``observable: <id>, formula: <expr>`` measurement models as
``{id: formula}`` (ADR-0036). A measurement model is a PEtab observableFormula evaluated
post-simulation by the observation layer; on export its ``.exp`` column classifies as the
measurement model (not a model entity), and its formula is emitted as the
``observableFormula`` verbatim -- the inverse of the importer's measurement-model line, so
an expression observable round-trips export -> import -> re-export."""
return {k[1]: v for k, v in conf.items()
if isinstance(k, tuple) and len(k) == 2 and k[0] == 'measurement'}
def _apply_observable_overrides(datas, overrides):
"""Rename each ``<header>`` data column (and its ``<header>_SD`` companion) to the
model ``<entity>`` across all ``datas`` so the column classifies against a model
observable/function -- mirroring ``config.py::_load_observables``. Global: a ``Data``
lacking a header is skipped (it just does not measure that observable); a header
present in **no** ``Data`` is a typo -> ``PybnfError`` (listing the columns present)."""
for entity, header in overrides.items():
found = False
for data in datas:
if header in data.cols:
data.rename_column(header, entity)
found = True
sd = f'{header}_SD'
if sd in data.cols:
data.rename_column(sd, f'{entity}_SD')
found = True
if not found:
present = sorted({c for data in datas for c in data.cols})
raise PybnfError(
f"Observable override 'observable: {entity}, column: {header}' names data "
f"column '{header}', but no experimental data file contains it (columns "
f"present: {present}). Check for a typo in the column name.")
def _read_conditions(conf, models, registry):
"""Read + validate the new-era ``condition:`` entries from the raw ``ploop`` dict.
Each ``('condition', name)`` entry is ``(model_ref_or_None, [(var, op, val_str), ...])``
(a named set of parameter perturbations -- a PyBNF Mutant = a PEtab Condition). Returns
``(conditions, condition_models)``: ``conditions`` is ``{condition_name: [(var, op,
float(val)), ...]}`` and ``condition_models`` is ``{condition_name: model_file}``, the one
model the condition belongs to.
A PEtab condition is model-agnostic (no modelId column; ADR-0041), but a PyBNF condition
belongs to exactly one model (ADR-0041 addendum): its ``model:`` ref, or the sole model of
a single-model job. The fitter attaches the condition to that model only and reads a fixed
target's base value from it (``bngsim_model/expressions.py::_nominal_param_value``), so the
exporter resolves the same model here -- a multi-model condition with no ``model:`` ref is
refused exactly as ``config.py::_load_conditions`` refuses it -- and validates each
perturbation target against **that** model's parameters / compartments. Two models may
give a fixed parameter of the same name different values, so the model a relative op is
computed against matters (#897). The single-model job validates against its one model
exactly as before."""
stem_to_model = {Path(mf).stem: mf for mf in models}
conditions = {}
condition_models = {}
for key, value in conf.items():
if not (isinstance(key, tuple) and len(key) == 2 and key[0] == 'condition'):
continue
name = key[1]
model_ref, perts = value
if model_ref is not None:
if Path(model_ref).stem not in stem_to_model:
raise PybnfError(
f"Condition '{name}' is declared for model '{model_ref}', but the job "
f"declares no model with id '{Path(model_ref).stem}' (declared model ids: "
f"{sorted(stem_to_model)}).")
model_file = stem_to_model[Path(model_ref).stem]
elif len(models) == 1:
model_file = models[0]
else:
raise PybnfError(
f"Condition '{name}' does not name a model, but the job declares {len(models)} "
f"models ({models}). Add 'model: <file>' to the condition to say which model it "
f"perturbs (ADR-0041); the fitter refuses this condition for the same reason.")
condition_models[name] = model_file
view = registry[model_file]
muts = []
for var, op, val in perts:
# A species-target perturbation (setConcentration -- a wash/bolus, #474) has a BNGL
# pattern target (contains '(') and possibly a param-expression value; its PEtab
# export is deferred, so pass it through UNVALIDATED (no param/compartment check, no
# float) -- the deferred-export guard in _build_experiments_and_conditions raises a
# clear message for a *referenced* one, and an unused species condition is skipped.
if '(' in var:
muts.append((var, op, val))
continue
if var not in view.parameters and var not in view.compartment_names:
where = (f"model '{model_file}', the model the condition belongs to"
if len(models) > 1 else 'the model')
raise PybnfError(
f"Condition '{name}' perturbs '{var}', which is not a parameter or "
f"compartment of {where} (a PEtab condition target must be a model "
f"entity).")
try:
muts.append((var, op, float(val)))
except (TypeError, ValueError):
# A parameter-valued perturbation (a per-condition estimated initial condition,
# ADR-0076): the value names a free parameter, passed through as a string; the
# builder emits it verbatim as the PEtab targetValue (mutation_target_value).
muts.append((var, op, val))
conditions[name] = muts
return conditions, condition_models
# ---------------------------------------------------------------------------
# Scope resolution (the documented boundaries)
# ---------------------------------------------------------------------------
def _resolve_models(conf):
"""Return the job's model files in declaration order (ADR-0041).
A BNGL (``.bngl``) or SBML (``.xml``) model is exported in its own native language
(ADR-0040, dispatched by :func:`_model_language`); any other extension raises there. A
job may declare one or many models; the new-era ``model:`` declarations accumulate into
``conf['model']`` in declaration order (a legacy ``model = X : Y`` job has none -- it is
refused downstream by :func:`_require_new_era_data` -- so fall back to the model set).
The model id is the file **stem** (the ``model_files`` key); two files sharing a stem
would collide on that key and the output filename, so a stem collision raises."""
ordered = list(dict.fromkeys(conf.get('model', [])))
if not ordered:
ordered = sorted(conf.get('models', set()))
if not ordered:
raise PybnfError(
"The job declares no model; add a 'model: <file>' declaration (ADR-0028).")
stems = {}
for mf in ordered:
_model_language(mf) # validate the extension is a supported language
stem = Path(mf).stem
if stem in stems:
raise PybnfError(
f"Models '{stems[stem]}' and '{mf}' share the file stem '{stem}', which "
f"would collide on the PEtab modelId and the exported filename. Rename one "
f"so each model has a distinct stem (ADR-0041).")
stems[stem] = mf
return ordered
def _require_modern_edition(conf):
"""Refuse a legacy-edition job: PEtab v2 interop is a new-era (``edition >= 2``)
feature (Bill's call, ADR-0031). A legacy conf names its objective with the retired
``objfunc`` key and binds data through the filename->suffix linkage; the exporter
reads only the modern surface, so it requires the conf to have opted into the new era
rather than reverse-mapping legacy syntax. Gates the *exporter* only -- the fitter
still runs legacy confs unchanged."""
ed = edition.resolve_edition(conf.get('edition'))
if not edition.is_modern(ed):
raise NotImplementedError(
"The PEtab v2 exporter requires a new-era config (edition >= 2); this job is "
f"legacy (edition {ed}). PEtab v2 interop is a new-era feature: add "
f"'edition = {edition.CURRENT_EDITION}' and name the objective on the modern "
"surface ('objective = <name>' or 'noise_model = <family>, ...') instead of "
"the legacy 'objfunc' key (ADR-0031, #423).")
def _resolve_noise(conf):
"""The job's whole-fit noise model as ``(noiseDistribution, sigma_verb, sigma_arg)``.
Modern-only (``export_job`` has already required ``edition >= 2``): the objective is
named on the ADR-0031 surface -- a whole-fit ``noise_model = <family>, ...`` line or
the named ``objective`` token -- with **no legacy** ``objfunc`` and **no implicit
default**, mirroring ``config.py``'s modern ``_load_obj_func`` branch. The resolved
objective is reduced to one ``(family, {param: (verb, arg)}, location)`` tuple and
reversed to PEtab.
The result is the **whole-fit base** applied to every column that has no per-observable
``noise_model <obs> = ...`` override; the overrides are resolved separately by
:func:`_resolve_per_observable_noise` (the additive seam, ADR-0021/0045) and a column
with one uses its own sigma source instead.
Raises ``NotImplementedError`` at every PEtab boundary, never a silent default:
* a column-joint ``profile_objective`` (``kl`` / ``wasserstein``) -- it scores the
whole column's shape, not a per-observation likelihood, so it has no PEtab
observable-noise representation;
* no objective, or more than one global objective key (no implicit default);
* a ``mean``-centered noise model on a family whose mean is not its median
(``lnnormal``) -- PEtab takes the prediction as the median for every family. The
location is the one the fitter uses: the global ``noise_location`` key when set (it
overrides the line's own ``location`` field, #898), else the line's field;
* an objective with no per-point noise model (``score`` / unknown token);
* a family PEtab v2 cannot express (``neg_bin`` -- removed; ``lognormal`` -- log10 vs
PEtab natural log). The distinct ``lnnormal`` family maps exactly to ``log-normal``.
"""
if conf.get('profile_objective') is not None:
raise NotImplementedError(
f"profile_objective = {conf['profile_objective']!r} is a column-joint "
f"objective (kl / wasserstein): it scores the whole column's shape, not a "
f"per-observation likelihood, so it has no PEtab observable-noise "
f"representation (ADR-0031, #423).")
whole_fit = conf.get(('noise_model', None))
has_objective = conf.get('objective') is not None
if whole_fit is not None and has_objective:
raise PybnfError(
"Specify exactly one global objective: this job has both a whole-fit "
"'noise_model = ...' line and an 'objective = ...' key.")
if whole_fit is not None:
family_token, fields, location = whole_fit # modern whole-fit line
elif has_objective:
token = conf['objective']
if token not in _OBJECTIVE_DESUGAR:
raise NotImplementedError(
f"objective = {token!r} has no per-point PEtab noise model: 'score' (no "
f"likelihood) and any unknown token are not PEtab observable noise. "
f"Per-point objectives: {sorted(_OBJECTIVE_DESUGAR)} (ADR-0031, #423).")
family_token, fields, location = _OBJECTIVE_DESUGAR[token](conf)
else:
raise NotImplementedError(
"No objective is named. A new-era (edition >= 2) job must name its objective "
"explicitly -- there is no implicit default. Set 'objective = <name>' or "
"'noise_model = <family>, ...' (ADR-0031, #423).")
# The global ``noise_location`` key (ADR-0024) is the whole-fit default location, and the
# fitter applies it last: ``Configuration._load_obj_func`` builds the objective from the
# line or token above, then calls ``set_default_location``, which rebuilds the class-default
# noise model with this location whatever the line's own ``location`` field said. It
# reaches only that class default -- the model every column without a per-observable
# ``noise_model <obs> = ...`` override is scored with -- never an override, which keeps its
# own location (``_resolve_per_observable_noise``). The exporter read only the line or
# token, so a mean-centred ``lnnormal`` fit was written as PEtab's median-centred
# ``log-normal`` with no refusal (#898). Reading the key here makes the location the
# exporter judges the one the fitter scores with.
where = 'the whole-fit noise model'
global_location = conf.get('noise_location')
if global_location is not None:
if global_location not in ('mean', 'median'):
raise PybnfError(
f"noise_location must be 'mean' or 'median', not {global_location!r}.")
location = global_location
where = f"the whole-fit noise model (noise_location = {global_location})"
return _reduce_noise_spec(family_token, fields, location, where)
def _resolve_per_observable_noise(conf):
"""The per-observable ``noise_model <obs> = ...`` overrides as ``{column:
(noiseDistribution, sigma_verb, sigma_arg)}`` (ADR-0021/0045) -- the additive companion to
:func:`_resolve_noise`'s whole-fit base.
Each ``('noise_model', <col>)`` key (``<col>`` not ``None``) is the parsed
``(family_token, {param: (verb, arg)}, location)`` spec for one observable **column** (the
model entity / measurement-model column the objective scores), reduced through the same
:func:`_reduce_noise_spec` boundaries as the whole-fit case (a ``mean`` location or a family
PEtab v2 cannot express raises). Empty when the job declares no override -- then every column
takes the whole-fit base and the export is byte-identical to the pre-per-observable output.
The inverse of the importer's :func:`~pybnf.petab.import_._per_observable_directives`, the
config side that consumes these (``objective._build_noise_overrides``)."""
overrides = {}
for key, value in conf.items():
if isinstance(key, tuple) and len(key) == 2 and key[0] == 'noise_model' \
and key[1] is not None:
column = key[1]
family_token, fields, location = value
overrides[column] = _reduce_noise_spec(
family_token, fields, location, f"the 'noise_model {column}' override")
return overrides
def _reject_cumulative(conf):
"""Fail loud if the job declares a ``cumulative`` prediction transform (ADR-0051, #418).
The cumulative->incident differencing is a PyBNF *prediction* transform with no PEtab v2
representation -- PEtab observables/measurements have no row-coupled cumulative-counts
operator. Exporting would silently drop it and emit a problem that scores the raw
cumulative columns, a different objective. Refuse instead (the project's fail-loud-over-
silently-wrong stance), naming the offending observables."""
cumulative = sorted(k[1] for k in conf
if isinstance(k, tuple) and k[0] == 'cumulative')
if cumulative:
raise NotImplementedError(
f"Observable(s) {cumulative} declare a cumulative->incident prediction transform "
f"('cumulative', #418), which PEtab v2 cannot express -- it has no row-coupled "
f"cumulative-counts observable operator. Exporting would silently score the raw "
f"cumulative columns instead. Remove the 'cumulative' flag (and difference the data "
f"to per-interval increments yourself) to export to PEtab.")
def _reject_time_error(conf):
"""Fail loud if the job marginalizes the latent measurement time (ADR-0112, #587).
A ``time_error`` clause on a ``noise_model`` line says the reported times are **not
exact**: the objective integrates each observation's density over a prior on its true
sampling time, and ``config.py`` swaps the whole per-point objective for a
``MarginalizedTimeObjective`` (a fit that reads ``MarginalizedTimeObjective`` before an
export reads ``LikelihoodObjective`` after one). A PEtab ``measurements`` row carries one
exact ``time`` and has no field for a distribution over it, so there is nothing to write
the clause as: exporting emits a problem scored at the nominal reported times by the
ordinary likelihood, which is a different statistical model.
Refuse instead -- the same stance and the same reason as :func:`_reject_cumulative`, whose
clause is a sibling in the very same ``noise_model`` grammar and was already guarded here
while this one was not (#738). The ``sigma_t`` scale source rides on the same key, so it
goes with it.
"""
keys = [k for k in conf if isinstance(k, tuple) and k[0] == 'time_error']
if not keys:
return
# The observable is None for the whole-fit form, which is the only one a job can
# currently run (``_maybe_marginalize_time`` defers per-observable time priors), but the
# exporter reads the raw config, so name whichever shape is actually present.
observables = sorted(k[1] for k in keys if k[1] is not None)
subject = (f"Observable(s) {observables} declare" if observables
else "This fit declares")
raise NotImplementedError(
f"{subject} a 'time_error' measurement-time marginalization (ADR-0112, #587), which "
f"PEtab v2 cannot express -- a measurements row carries one exact 'time' and has no "
f"field for a distribution over it. Exporting would emit a problem scored at the "
f"nominal reported times by the ordinary likelihood, a different objective. Remove "
f"the 'time_error' and 'sigma_t' fields from the noise_model line to export to "
f"PEtab.")
def _reject_normalization(conf):
"""Fail loud if the job declares any normalization (ADR-0053, #444; floor/scale ADR-0066, #479).
Normalization (``peak`` / ``init`` / ``zero`` / ``unit`` / ``floor`` / ``scale``) is a PyBNF
*prediction* transform -- a whole-trajectory reduction of a predicted observable before scoring
-- with no PEtab v2 representation: PEtab observable formulas are pointwise, so they cannot
express "divide by this trajectory's peak / initial value" (the scale comes from the trajectory
itself, not the model state at one point). ``floor`` (``x + rho*max(x)``) is a whole-series
offset and ``scale`` is an analytic per-series optimum profiled at scoring time -- both equally
non-pointwise (``scale`` does have a natural PEtab home in estimated ``observableParameters`` /
hierarchical scaling, a future export mapping, #479). Exporting would silently drop the
transform and emit a problem that scores the raw, un-normalized columns, a different objective.
Refuse instead (the fail-loud-over-silently-wrong stance, like :func:`_reject_cumulative`),
naming what is normalized -- the per-observable ``('normalization', target)`` keys (ADR-0053),
the whole-fit / legacy ``normalization`` value, and the compiled ``analytic_scale`` key (which
a whole-fit ``scale`` reduces to with no ``normalization`` value left)."""
targets = sorted(k[1] for k in conf
if isinstance(k, tuple) and k[0] == 'normalization')
whole_fit = conf.get('normalization')
analytic_scale = conf.get('analytic_scale')
if not targets and whole_fit is None and not analytic_scale:
return
if targets:
detail = f"observable(s) {targets}"
elif whole_fit is None and analytic_scale:
cols = sorted({c for cols in analytic_scale.values() for c in cols})
detail = f"observable(s) {cols} (analytic per-series scaling)"
elif isinstance(whole_fit, dict):
detail = f"data file(s) {sorted(whole_fit)}"
else:
detail = f"the whole fit ('{whole_fit}')"
raise NotImplementedError(
f"This job normalizes {detail} (the 'normalization' key, ADR-0053/ADR-0066), a PyBNF "
f"prediction transform PEtab v2 cannot express -- it has no observable operator for "
f"peak/initial-value/z-score/floor normalization or analytic per-series scaling (a "
f"whole-trajectory reduction, not a pointwise observable formula). Exporting would "
f"silently score the raw, un-normalized columns instead. Remove the normalization "
f"(normalizing your data and model output equivalently yourself) to export to PEtab.")
def _reject_postprocess(conf):
"""Fail loud if the job runs a ``postprocess`` script on its simulations (#899).
``postprocess = <script.py> <suffix> ...`` names a user Python function that the fitter
applies to each named simulation before scoring it (``Configuration._load_postprocessing``
maps the ``(model, suffix)`` pairs, under edition 2 the experiment name is the suffix;
``Result.postprocess_data`` swaps the simulation for ``postprocess(data)``). An arbitrary
Python transform of the prediction has no PEtab v2 representation -- an ``observableFormula``
is a pointwise expression over model entities, not a program -- so exporting wrote the bare
model column and emitted a problem that scores the untransformed simulation: a different
objective with a different optimum, and no warning (#899). Refuse instead, beside
:func:`_reject_normalization`, whose whole-trajectory reductions are the built-in cousins of
what a script typically does. The key is model-language agnostic (the fitter applies it to
a BNGL or an SBML model's suffix alike), so is the refusal."""
specs = conf.get('postprocess')
if not specs:
return
detail = '; '.join(
f"'{spec[0]}' on {', '.join(repr(s) for s in spec[1:])}" for spec in specs)
raise NotImplementedError(
f"This job transforms its simulations with a 'postprocess' script ({detail}) before "
f"scoring them. A user Python function applied to the prediction has no PEtab v2 "
f"representation -- an observableFormula is a pointwise expression over model "
f"entities, not a program -- so exporting would silently score the untransformed "
f"simulations instead, a different objective (#899). Remove the 'postprocess' line to "
f"export, after expressing the transform as an 'observable: <id>, formula: <expr>' "
f"measurement model if it is pointwise, or applying it to the data yourself.")
def _reduce_noise_spec(family_token, fields, location, where):
"""Reduce one parsed noise spec ``(family_token, {param: (verb, arg)}, location)`` to
``(noiseDistribution, sigma_verb, sigma_arg)``, raising the PEtab boundaries shared by the
whole-fit base (:func:`_resolve_noise`) and the per-observable overrides
(:func:`_resolve_per_observable_noise`): a ``mean``-centered location on a family whose
mean is not its median (PEtab is median-only) and a family PEtab v2 cannot express
(``neg_bin`` removed; ``lognormal`` is log10 vs PEtab's natural ``log-normal``; ``lnnormal``
is its exact native match). ``where`` names the spec in the error message.
A ``mean`` location on a linear Gaussian or Laplace exports: those families are symmetric,
so the mean IS the median and the fitter's likelihood is PEtab's to the bit (the moment
offset is exactly ``0.0``; :data:`_MEAN_IS_MEDIAN_FAMILIES`). Refusing it would only make
the user delete a word that changes no number (#898)."""
if location == 'mean' and family_token.lower() not in _MEAN_IS_MEDIAN_FAMILIES:
raise NotImplementedError(
f"{where} is mean-centered (location = mean) on the '{family_token}' family, "
f"whose mean is not its median; PEtab v2 takes the prediction as the distribution "
f"median for every noise family, so mean centering has no PEtab representation "
f"and the exported problem would have a different optimum (ADR-0031, #423, #898). "
f"Use median (drop 'noise_location = mean' or the line's 'location = mean').")
distribution = _FAMILY_TOKEN_TO_PETAB_DISTRIBUTION.get(family_token.lower())
if distribution is None:
raise NotImplementedError(
f"{where}: the '{family_token}' noise family cannot be expressed in PEtab v2: "
f"neg_bin was removed from v2, and PyBNF's lognormal is log10 while PEtab's "
f"log-normal is natural log (use the distinct lnnormal family for that scale). "
f"ADR-0023/0031, #423.")
(_param, (verb, arg)), = fields.items()
return distribution, verb, arg
def _independent_variable(data):
"""The header of a wide :class:`~pybnf.data.Data`'s column 0 (``time`` or the swept axis)."""
return min(data.cols, key=data.cols.get)
def _swept_param(exp):
"""The swept-parameter (dose axis) header of a parameter_scan experiment's data (column 0)."""
data0 = exp['datas'][0]
return data0.indvar if data0.indvar is not None else _independent_variable(data0)
def _refuse_wash_re_pins(pdr_experiments, conditions, surrogate):
"""Refuse a pre-equilibrated dose-response whose wash condition would re-pin a parameter it
must leave alone (#892).
A wash (measurement) condition is emitted like every other condition: it carries a base pin
``p = p__REF`` for each ``p`` in the surrogate set M that it does not set itself (ADR-0027).
On this shape two such pins are wrong:
* ``p`` is set by the pre-equilibration condition. The fitter applies that condition as an
inline ``setParameter`` and never undoes it, so ``p`` keeps the pre-equilibration value
through the dose scan. PEtab v2 keeps it too, unless a later period sets it again, which is
exactly what the wash's pin would do.
* ``p`` is the swept parameter. The per-dose condition sets it in the same period, and PEtab
v2 forbids two conditions of one period from setting the same target.
For the same reason a wash that sets the swept parameter itself is refused, whether or not
the parameter is in M. The fitter's scan overrides the wash's value at every dose, but the
exported measurement period would give the swept parameter two setters.
A wash-free scan has no such pin (its measurement period carries only the per-dose
condition), so it is never refused here.
"""
for exp in pdr_experiments:
wash = exp['condition']
if wash is None:
continue
wash_targets = {var for var, _op, _val in conditions[wash]}
swept = _swept_param(exp)
if swept in wash_targets:
raise NotImplementedError(
f"Pre-equilibrated dose-response experiment '{exp['name']}' sweeps '{swept}', and "
f"its wash condition '{wash}' also sets '{swept}'. Both conditions apply in the "
f"exported measurement period, and PEtab v2 forbids two conditions of one period "
f"from setting the same target. PyBNF's scan sets '{swept}' to each dose over the "
f"wash's value, so remove '{swept}' from '{wash}' to export this job, or run it "
f"natively.")
pre = exp['preequilibrate']
# The swept parameter is left to the check below: the dose, not the pre-equilibration
# value, is what the scan runs at.
carried = sorted(p for p in surrogate
if p != swept and p not in wash_targets
and any(var == p for var, _op, _val in conditions[pre]))
if carried:
raise NotImplementedError(
f"Pre-equilibrated dose-response experiment '{exp['name']}': its "
f"pre-equilibration condition '{pre}' sets the fit parameter(s) {carried}, and "
f"PyBNF keeps that value through the dose scan (the pre-equilibration setParameter "
f"is never undone). In the exported problem the wash condition '{wash}' must re-pin "
f"every fit-and-perturbed parameter it does not set to its estimate "
f"({carried[0]} = {surrogate_name(carried[0])}, ADR-0027), which would undo that "
f"value. Repeat the pre-equilibration value in '{wash}' so the scan's value is "
f"stated there too, or run the job natively.")
if swept in surrogate:
raise NotImplementedError(
f"Pre-equilibrated dose-response experiment '{exp['name']}' sweeps '{swept}', a fit "
f"parameter that a condition also perturbs, so its wash condition '{wash}' re-pins "
f"'{swept}' to its estimate ({swept} = {surrogate_name(swept)}, ADR-0027) in the "
f"same period in which the per-dose condition sets it to the dose. PEtab v2 forbids "
f"two conditions of one period from setting the same target. Fix '{swept}' (do not "
f"fit it) to export this job, or run it natively.")
# The relative tolerance the fitter matches a data row's dose with: ``Objective._sim_row_for``
# takes the first simulated row where ``np.isclose(sim_dose, dose, atol=0.)`` holds, i.e. numpy's
# default ``rtol``. Two distinct doses this close can be scored against one simulation.
_FITTER_DOSE_RTOL = 1e-5
def _dose_axis(exp):
"""The exported dose axis of a parameter_scan experiment: the sorted union of every
replicate's doses, as floats (#895).
This is exactly the grid the fitter scans (``config.py`` stacks the ``data:`` replicates and
passes ``sorted({float(x) for x in stacked[indvar]})`` to the scan), so the exported problem
has one Experiment per dose the fitter simulates. The exporter then tags each data row with
the experiment of its own dose (:func:`~pybnf.petab.measurements.dose_response_measurement_rows`),
as the fitter pairs each row with the simulation at its own dose. The swept column is read by
name, as the fitter's replicate stacking reads it.
Refuses (``PybnfError``, naming the experiment, file and dose) a data file the fitter could not
score correctly either: a replicate with no swept-parameter column (its doses stack as NaN); a
non-finite dose; and two distinct doses within the fitter's dose-matching tolerance
(:data:`_FITTER_DOSE_RTOL`), which the fitter may score against a single simulation while PEtab
would simulate each separately.
"""
swept = _swept_param(exp)
first_file = {}
for data, data_file in zip(exp['datas'], exp['data_files']):
if swept not in data.cols:
raise PybnfError(
f"Experiment '{exp['name']}' is a dose-response scan over '{swept}' (the first "
f"column of its first data file), but data file '{data_file}' has no '{swept}' "
f"column (its columns: {list(data.cols)}). Every replicate of a dose-response "
f"must give each row's dose in a '{swept}' column.")
for value in data[swept]:
dose = float(value)
if not np.isfinite(dose):
raise PybnfError(
f"Experiment '{exp['name']}', data file '{data_file}': the dose "
f"{swept} = {dose!r} is not a finite number, so no simulation can be run or "
f"matched at it. Give every row a finite dose.")
first_file.setdefault(dose, data_file)
axis = sorted(first_file)
for low, high in zip(axis, axis[1:]):
if np.isclose(low, high, rtol=_FITTER_DOSE_RTOL, atol=0.0):
raise PybnfError(
f"Experiment '{exp['name']}' has two distinct doses {swept} = {low!r} (in "
f"'{first_file[low]}') and {swept} = {high!r} (in '{first_file[high]}') that "
f"differ by less than the fitter's dose-matching tolerance (relative "
f"{_FITTER_DOSE_RTOL:g}). The fitter can score measurements at one of them against "
f"the simulation of the other, while the exported PEtab problem would simulate "
f"each separately. If they are the same dose, write it identically in every data "
f"file.")
return axis
def _species_id_map(patterns):
"""A ``{pattern: petab_id}`` map for the species ``setConcentration`` targets the job's
conditions reference (ADR-0062), in first-appearance order. Each BNGL pattern is aliased to a
synthesized PEtab id (:func:`~pybnf.petab.conditions.species_target_id`) for the mapping table;
two distinct patterns that sanitize to the same id collide -> ``PybnfError``."""
mapping = {}
by_id = {}
for pattern in patterns:
if pattern in mapping:
continue
pid = species_target_id(pattern)
if pid in by_id:
raise PybnfError(
f"Species setConcentration targets {by_id[pid]!r} and {pattern!r} both sanitize "
f"to the PEtab mapping id {pid!r}. Rename one species pattern so their PEtab "
f"target ids stay distinct (ADR-0062).")
by_id[pid] = pattern
mapping[pattern] = pid
return mapping
# ---------------------------------------------------------------------------
# Observable + parameter rows
# ---------------------------------------------------------------------------
def _observable_rows(experiments, registry, noise, per_obs_noise, inline_functions=False,
measurement_models=None):
"""Classify each fitted column across all experiments as a model observable, a model
function, or a conf measurement model, and map it to a PEtab observable row.
A column is classified against the model of **each** experiment that measures it
(ADR-0041): the observables table has no per-model namespace, so a column shared across
models must classify identically in each (the same kind, hence the same observableId and
formula) -- a column that is an observable in one model and a function in another is a
real conflict and raises. A column is gathered once, in first-appearance order across the
experiments' (override-renamed) ``datas``, so the observables table covers the whole job.
``noise`` is the whole-fit base ``(noiseDistribution, sigma_verb, sigma_arg)`` from
:func:`_resolve_noise` and ``per_obs_noise`` the ``{column: (dist, verb, arg)}``
per-observable overrides (ADR-0021/0045): each column's noise (its family + sigma source)
is its override if present, else the base, resolved across every experiment's data (it can
depend on the column's data, e.g. a ``column_mean`` sigma).
``inline_functions`` (ADR-0035) emits a **function** column's body as an
``observableFormula`` expression instead of the bare name -- the opt-in path that
generates the importer's round-trip oracle; the default keeps every column bare.
``measurement_models`` (``{id: formula}``, ADR-0036) are conf-declared measurement models
(model-agnostic): a column matching one is emitted with that formula as its
``observableFormula`` and its id verbatim (the inverse of the importer's ``observable: ...
formula:`` line).
Returns ``(observable_rows, column_to_observable_id, column_means)``. ``column_means`` is
``{column: {experiment name: mean}}`` for each ``column_mean`` column whose experiments
have different means (#894): its observable row declares a noise placeholder, and the
caller writes each measurement row's own experiment mean into ``noiseParameters``."""
measurement_models = measurement_models or {}
# Gather the fitted columns in first-appearance order, each tagged with the model
# file(s) that measure it (distinct, declaration order). A measurement-model column is
# model-agnostic (it never touches a model's entities), so its models list is inert.
columns = []
column_models = {}
for exp in experiments:
mf = exp['model']
for data in exp['datas']:
indvar = data.indvar if data.indvar is not None else _independent_variable(data)
for col in sorted(data.cols, key=data.cols.get):
if col == indvar or col.endswith('_SD'):
continue
if col not in column_models:
columns.append(col)
column_models[col] = []
if mf not in column_models[col]:
column_models[col].append(mf)
observable_rows = []
column_to_observable_id = {}
column_means = {}
for col in columns:
classes = [(mf, _classify_column(col, registry[mf], mf, measurement_models,
inline_functions))
for mf in column_models[col]]
distinct = {kind_formula for _mf, kind_formula in classes}
if len(distinct) != 1:
detail = '; '.join(f"{mf} -> {kind}" for mf, (kind, _f) in classes)
raise PybnfError(
f"Exp column '{col}' classifies inconsistently across the models that "
f"measure it ({detail}). PEtab's observables table has no per-model "
f"namespace, so a column shared across models must mean the same observable "
f"in each (ADR-0041).")
kind, formula = classes[0][1]
# A column's noise is its per-observable override if one is declared, else the
# whole-fit base (ADR-0021/0045); the override carries its own family + sigma source.
distribution, verb, arg = per_obs_noise.get(col, noise)
noise_source = _noise_source_for_column(verb, arg, col, experiments)
if noise_source[0] == 'experiment_means':
# A column-mean sigma that differs between experiments (#894): the observable
# declares a noise placeholder and each measurement row carries its own
# experiment's mean in noiseParameters -- the fit's per-experiment sigma, exactly.
column_means[col] = noise_source[1]
noise_source = ('placeholder', None)
row = petab_observable_row(col, kind, distribution, noise_source,
observable_formula=formula)
observable_rows.append(row)
column_to_observable_id[col] = row.observable_id
if not observable_rows:
raise PybnfError(
"The job's experiment data has no fittable observable/function columns "
"(only an independent variable and/or _SD columns).")
return observable_rows, column_to_observable_id, column_means
def _classify_column(col, model, model_file, measurement_models, inline_functions):
"""Classify one fitted column against one model view, returning ``(kind, formula)``
(the formula is ``None`` unless it is a measurement model or an inlined function body).
A conf-declared measurement model wins first (it is model-agnostic, ADR-0036); else the
column must be a model observable or function (its bare name being the
``observableFormula``). A column matching nothing in this model raises -- naming the
model so a multi-model job points at the right namespace (ADR-0041)."""
if col in measurement_models:
return 'measurement', measurement_models[col] # the conf observableFormula, verbatim
if col in model.observable_names:
return 'observable', None
if col in model.function_names:
formula = _inlined_formula(col, 'function', model, model_file) \
if inline_functions else None
return 'function', formula
raise PybnfError(
f"Exp column '{col}' matches no observable, function, or measurement model in model "
f"'{model_file}' (its observables: {sorted(model.observable_names)}; functions: "
f"{sorted(model.function_names)}; measurement models: {sorted(measurement_models)}).")
def _inlined_formula(col, kind, model, model_file):
"""The ``observableFormula`` for a column under inlining mode (ADR-0035), or ``None``.
Only a **function** column is inlined (an observable is a model species/group, not an
algebraic expression, so it stays bare); its captured body is translated to PEtab math
by :func:`~pybnf.petab.formula.bngl_body_to_petab_math`. A function with an empty body
-- a forward declaration, or a function *of arguments* (only zero-arg global functions
are the BNGL measurement-model convention) -- cannot be inlined and raises rather than
emitting a bare-name formula that silently contradicts the requested mode. Reached only
for a BNGL model (SBML has no functions -- :class:`_SbmlModelView.function_names` is
empty), so ``model`` here is always a :class:`~pybnf.petab._bngl.BnglEntities`.
"""
if kind != 'function':
return None
body = model.function_bodies.get(col, '')
if not body:
raise NotImplementedError(
f"Function '{col}' in model '{model_file}' has no inlinable body (a forward "
f"declaration or a function with arguments); only a zero-arg global function "
f"'{col}() = <body>' can be inlined as an observableFormula (ADR-0035). Export "
f"without inline_functions to reference it by bare name.")
# The body is BNGL, read with BioNetGen's grammar, not PEtab's (#908).
return bngl_body_to_petab_math(body, model, function_name=col, model_file=model_file)
def _noise_source_for_column(verb, arg, col, experiments):
"""The PEtab noise representation for one fitted column, from the desugared sigma
source verb (ADR-0021 reversed) -- evaluated across every experiment's ``datas``:
* ``read_exp_file`` (the ``_SD`` data column) -> a per-point placeholder, fed by the
measurements' ``noiseParameters``. Every ``Data`` carrying the column must also
carry its ``<col><suffix>`` companion (else a measurement row would lack the noise
value its declared placeholder binds to).
* ``fix_at`` -> a constant noiseFormula (the fixed sigma).
* ``column_mean`` -> the column's mean **per experiment** (:func:`_column_mean_noise_source`,
#894): a constant noiseFormula when every experiment measuring the column has the same
mean (always so for one experiment), else ``('experiment_means', {name: mean})``, which
the caller turns into a noise placeholder fed row by row.
* ``formula`` -> the expression noiseFormula verbatim (a ``FormulaSigma``, ADR-0044/0045):
a PEtab-math expression over free-parameter ids + constants. The expression's symbols are
PEtab parameter ids (exported as estimated parameters); a noise nuisance that is not a
model parameter is still a deferred boundary (it would fail the model-id binding check,
shared with the ``fit`` sigma), so the whole-fit ``formula`` export covers an expression
over model parameters. When the ``formula`` expression still carries a per-measurement
**placeholder** (``noiseParameter*``), the sigma is row-varying
(``PerMeasurementFormulaSigma``, ADR-0045): it returns ``('per_measurement', expr)``, the
noiseFormula emitted verbatim with its placeholder and the per-row token supplied by the
measurements' ``noiseParameters`` column (the binding-table sidecar).
* ``fit`` -> ``('free_param', id)``: a free-parameter (estimated) sigma -> a bare-id
noiseFormula naming the noise parameter ``id``, declared estimated in the parameter table
and admitted as an observation-layer nuisance (not a model entity -- #439). The importer
reads it back to a ``fit`` source (ADR-0044), so a per-observable estimated sigma
round-trips.
A relative sigma (``relative``) is still a deferred boundary: it is a ``noiseFormula``
expression (the sympy layer, mirroring the importer's expression boundary).
"""
holders = [data for exp in experiments for data in exp['datas'] if col in data.cols]
if verb == 'formula':
return ('per_measurement', arg) if _PLACEHOLDER.search(arg) else ('formula', arg)
if verb == 'prediction_formula':
# A prediction-dependent sigma (PredictionFormulaSigma, ADR-0075): sigma scales with the
# simulated output, e.g. the combined error model ``sd_abs + sd_rel*y``. It maps back to
# a plain noiseFormula emitted VERBATIM -- the direct mirror of ``formula`` -- because the
# importer reclassifies it: the same substituted expression re-imports as
# ``prediction_formula`` iff it references a model entity (``y``), else ``formula``
# (_resolve_noise, import_.py). So the exporter carries no separate prediction arm: the
# coefficients (``sd_abs`` / ``sd_rel``) are admitted as nuisances, and the model-entity
# symbol stays a model reference in the noiseFormula (fit-preserving round trip, #502).
return ('formula', arg)
if verb == 'read_exp_file':
sd_col = col + arg
if any(sd_col not in data.cols for data in holders):
raise NotImplementedError(
f"Observable column '{col}': the objective reads its noise from the "
f"'{sd_col}' data column, but a data file carrying '{col}' has no such "
f"column. A constant or free-parameter sigma without per-point data is a "
f"separate path (ADR-0023, #423).")
return ('placeholder', None)
if verb == 'fix_at':
return ('constant', float(arg))
if verb == 'column_mean':
return _column_mean_noise_source(col, experiments)
if verb == 'fit':
# A free-parameter (estimated) sigma -> a bare-id noiseFormula naming the noise
# parameter (declared estimated in parameters.tsv; admitted as an observation-layer
# nuisance by _referenced_nuisance_symbols, NOT a model entity -- #439). The importer
# reads a bare-id noiseFormula back to a 'fit' source (ADR-0044), so a per-observable
# estimated sigma (the per_observable_noise example / Boehm's sd_*) round-trips.
return ('free_param', arg)
raise NotImplementedError(
f"Observable column '{col}': the '{verb}' sigma source is a later export chunk "
f"-- a relative sigma is a noiseFormula expression (the sympy layer, mirroring "
f"the importer boundary). ADR-0021/0023, #423.")
def _experiment_column_means(col, experiments):
"""``{experiment name: mean}`` -- the ``column_mean`` sigma the fit gives column ``col`` in
each experiment that has at least one observed value of it (#894).
This is the fit's own number. ``Objective.evaluate_multiple`` scores one experiment at a
time, and ``ColumnMeanSigma`` (like the legacy ``ave_norm_sos``) takes
``Data.column_mean`` of *that* experiment's Data. The Data is the experiment's replicate
files stacked in ``data:`` order (``config._stack_replicates``, ADR-0039), so its observed
values are these ``datas`` concatenated in the same order. The mean is therefore the same
float. It is taken over observed values only, since NaN means unmeasured (#707). A
dose-response is one PyBNF experiment even though its doses become N PEtab experiments, so
all of its doses share one mean. An experiment with no observed value of the column
contributes no scored point and no measurement row, so it has no sigma to export."""
means = {}
for exp in experiments:
values = [d[col] for d in exp['datas'] if col in d.cols]
if not values:
continue
mean = float(observed_mean(np.concatenate(values)))
if not np.isnan(mean):
means[exp['name']] = mean
return means
def _column_mean_noise_source(col, experiments):
"""The PEtab noise source for a ``column_mean`` sigma on column ``col`` (#894).
The fit normalizes each experiment by **its own** column mean
(:func:`_experiment_column_means`), so a single mean pooled over every experiment is the
wrong sigma as soon as two experiments measure the column at different magnitudes: it
reweights the experiments against each other and moves the optimum. PEtab has no
data-derived sigma, but the mean is a constant of the data, so it can be written exactly:
* every experiment has the same mean (in particular, only one experiment measures the
column) -> ``('constant', mean)``, the inline numeric noiseFormula (unchanged from the
pre-#894 export for a single experiment);
* the means differ -> ``('experiment_means', {name: mean})``: the observable declares a
noise placeholder (``noiseParameter1_<id>``) and every measurement row carries its own
experiment's mean in ``noiseParameters``. PEtab then gives each point the fit's sigma.
The importer reads either form back to ``column_mean`` only when each value equals its
experiment's own mean (``import_._ColumnMeans``), so the round trip restores the job.
"""
means = _experiment_column_means(col, experiments)
distinct = set(means.values())
if len(distinct) > 1:
return ('experiment_means', means)
# No observed value anywhere: no measurement row scores this column, so the constant is
# inert; keep the historical NaN (observed_mean of nothing) rather than invent a sigma.
return ('constant', distinct.pop() if distinct else float('nan'))
def _resolve_free_to_model(free_params, registry, models, nuisances=()):
"""Validate each free parameter binds to a model parameter id (or is an admitted
observation-layer nuisance); return the identity map.
New-era binds free parameters **by id** (ADR-0034): a free parameter's name *is* the
model parameter it drives -- no ``__FREE`` marker. The id space is the **union** of every
model's parameter ids (a BNGL ``begin parameters`` id, or an SBML global parameter id --
ADR-0040/0041); a free parameter present in at least one model binds (the same id means
the same knob, so a multi-model job binds it across all models that carry it). This is the
exporter's analogue of the new-era config typo check
(:meth:`config._check_variable_correspondence_modern`, which unions the same way): a free
parameter matching no model parameter id is a typo (or a ``fit`` sigma, which the exporter
rejects separately at column classification).
``nuisances`` are free parameters that bind to no model id but are legitimate
**observation-layer** nuisances (a measurement scale, a noise coefficient, a row-varying
per-row sigma -- ADR-0044/0045), gathered by :func:`_referenced_nuisance_symbols` from the
measurement-model / noise formulae and the binding-table tokens. They pass through as
estimated parameters with no model binding (the export peer of config widening the
measurement-model namespace + the ``_per_measurement_free_params`` orphan union); only a
free parameter that is neither a model id nor a referenced nuisance is a typo.
The exporter never builds a ``Configuration``, so it validates against the views'
``parameters`` directly. Returns ``{name: name}`` -- the identity map the rest of the
exporter threads as ``free_to_model``.
"""
model_ids = set().union(*(set(v.parameters) for v in registry.values()))
free_to_model = {}
for fp in free_params:
if fp.name not in model_ids and fp.name not in nuisances:
legacy_hint = ''
if fp.name.endswith('__FREE'):
legacy_hint = (
f" The '__FREE' marker is legacy-edition only (ADR-0034); declare "
f"the bare parameter id '{fp.name[:-len('__FREE')]}' instead.")
where = (f"model '{models[0]}'" if len(models) == 1
else f"any of the job's {len(models)} models ({models})")
raise PybnfError(
f"Free parameter '{fp.name}' matches no parameter id in {where}.",
f"Under edition >= 2 a BNGL free parameter binds to a model parameter by "
f"id (the SBML/PEtab convention; ADR-0034/0041), so '{fp.name}' must be one "
f"of the models' parameter ids: {sorted(model_ids)}.{legacy_hint}")
free_to_model[fp.name] = fp.name
return free_to_model
# A bare identifier in a measurement-model / noise formula: scanned to find which free
# parameters a formula references (over-matches model entities + placeholders, harmlessly --
# only the intersection with declared free parameters is used). Dependency-free (no petab).
_FORMULA_SYMBOL = re.compile(r'[A-Za-z_]\w*')
def _referenced_nuisance_symbols(conf, conf_path, noise, per_obs_noise):
"""The names referenced as observation-layer nuisances by the job's measurement-model and
noise surfaces (ADR-0034/0044/0045) -- the candidates :func:`_resolve_free_to_model` admits
as model-unbound estimated parameters. The union of:
* the symbols of every ``observable: <id>, formula: <expr>`` measurement-model formula
(an ``observableParameters`` scale/offset substituted in -- ADR-0044 -- reads as a free
symbol, e.g. ``scaling`` in ``scaling*x``);
* the symbols of every ``formula``-verb ``noiseFormula`` (whole-fit or per-observable),
e.g. ``slope`` in ``0.05*slope + 0.1`` (a ``FormulaSigma`` noise coefficient);
* the id of every ``fit``-verb sigma (whole-fit or per-observable), e.g. ``b_y`` in
``noise_model y = laplace, scale = fit b_y`` (a ``FreeParameterSigma`` estimated scale --
#439): an estimated noise parameter is an observation-layer nuisance, not a model entity;
* the non-numeric (parameter-id) tokens of every experiment's ``measurement_params:``
binding-table sidecar, e.g. the per-row ``sd_lo`` / ``s_lo`` (a row-varying estimated
sigma / scale -- ADR-0045).
Over-matches model entities and placeholders, harmlessly: the caller intersects with the
declared free parameters, so only a genuine free-parameter nuisance is admitted (an
unreferenced free parameter that is not a model id stays a typo)."""
referenced = set()
for formula in _read_measurement_models(conf).values():
referenced |= set(_FORMULA_SYMBOL.findall(formula))
for _dist, verb, arg in [noise, *per_obs_noise.values()]:
if verb in ('formula', 'prediction_formula'):
# A FormulaSigma expression (ADR-0044) or a prediction-dependent sigma (ADR-0075,
# ``sd_abs + sd_rel*y``): its free-parameter coefficients are nuisances. Over-matches
# the model-entity symbol (``y``) harmlessly -- the caller intersects with declared
# free parameters, which a model entity is not, so only the coefficients are admitted.
referenced |= set(_FORMULA_SYMBOL.findall(arg))
elif verb == 'fit':
referenced.add(arg) # an estimated noise scale (FreeParameterSigma), #439
for key, fields in conf.items():
if not (isinstance(key, tuple) and len(key) == 2 and key[0] == 'experiment'):
continue
mp_file = fields.get('measurement_params')
if not mp_file:
continue
table = read_measurement_params(conf_path.parent / mp_file)
for by_placeholder in table.values():
for by_time in by_placeholder.values():
referenced |= {tok for tok in by_time.values() if not _is_numeric_token(tok)}
# A free parameter a condition perturbation references by value (a per-condition estimated
# initial condition, ADR-0076): it binds no model entity of its own, so -- like a noise /
# measurement-model nuisance -- it must be admitted as a model-unbound estimated parameter.
for key, value in conf.items():
if not (isinstance(key, tuple) and len(key) == 2 and key[0] == 'condition'):
continue
_model_ref, perts = value
for var, _op, val in perts:
if '(' not in var and isinstance(val, str) and not _is_numeric_token(val):
referenced.add(val)
return referenced
def _is_numeric_token(token):
"""Whether a binding-table token is a numeric literal (inlined) vs a parameter id (an
estimated nuisance -- the export peer of ``config._is_numeric_token``)."""
try:
float(token)
return True
except (TypeError, ValueError):
return False
def _parameter_rows(free_params, free_to_model, surrogate_params, registry, models):
"""Map each free parameter to a row; a fit-and-mutated one renamed to ``<p>__REF``."""
union_ids = set().union(*(set(v.parameters) for v in registry.values()))
parameter_rows = []
for fp in free_params:
model_param = free_to_model[fp.name]
if model_param in surrogate_params:
ref = surrogate_name(model_param)
if ref in union_ids:
raise PybnfError(
f"The surrogate-base name '{ref}' for fit-and-mutated parameter "
f"'{model_param}' clashes with an existing parameter in the job's "
f"models. Rename that model parameter.")
parameter_id = ref
else:
parameter_id = model_param
parameter_rows.append(petab_parameter_row(fp, parameter_id=parameter_id))
return parameter_rows
def _condition_nominal_of(registry, condition_models):
"""The ``nominal_of(condition, var)`` callable the condition builders take: a fixed
parameter's numeric nominal value in the model **the condition belongs to**, or ``None``.
A relative op (``* / + -``) on a fixed target is folded to a number on export
(:func:`~pybnf.petab.conditions.mutation_target_value`), so the base it is folded against
must be the one the fitter uses: the condition's own model's value
(``bngsim_model/expressions.py::_nominal_param_value`` reads the experiment's engine model,
and an experiment can only apply a condition of its own model). Two models of a multi-model
job may give a same-named fixed parameter different values, so reading "the first model
that declares it" computed the condition against the wrong model (#897). A free target
never reaches here (the surrogate path handles it)."""
def nominal_of(condition, var):
return _numeric_nominal(registry[condition_models[condition]], var)
return nominal_of
def _numeric_nominal(model, var):
"""A fixed parameter's numeric nominal value, or ``None`` (expression/unknown RHS).
Works for both model views: a BNGL ``parameters`` value is the raw RHS string (a number
floats, an expression raises ``ValueError`` -> ``None``); an SBML view's value is already
a float or ``None`` (a value-less parameter -> ``None`` via the ``TypeError`` guard)."""
rhs = model.parameters.get(var)
if rhs is None:
return None
try:
return float(rhs)
except (ValueError, TypeError):
return None
# ---------------------------------------------------------------------------
# Reading the job (the disposable input half of the seam)
# ---------------------------------------------------------------------------
def _read_conf_dict(conf_path):
"""Parse a ``.conf`` to the raw ``ploop`` dict (no model loading, no BNG)."""
with open(conf_path) as fh:
return ploop(fh.readlines())
def _free_parameters_from_conf(conf):
"""Build ``FreeParameter`` objects from the config's free-parameter declarations.
Reads **both** spellings, in declaration order (ADR-0043): the legacy positional
``<family>_var = <id> p1 [p2]`` line and the new-era ``parameter:`` record. Only the
first was read until #733, so an edition-2 record was skipped rather than refused and
the whole free parameter vanished from the exported problem with no diagnostic --
including the truncated priors the *importer* emits as records, the one grammar
carrying ``lower``/``upper``.
Each parameter carries the fit's declared start point, if it has one, on ``.value`` --
the source :func:`~pybnf.petab.parameters.petab_parameter_row` writes as the row's
``nominalValue`` (#719). Both spellings of a start point are honoured here, the way
``Configuration._load_start_point`` merges them: a record's ``initial_value:`` field
and a ``start_point =`` line beside it.
"""
start_points = _start_points_from_conf(conf)
free_params = []
for key, value in conf.items():
if not (isinstance(key, tuple) and len(key) == 2
and isinstance(key[0], str) and isinstance(key[1], str)):
continue
keyword, name = key
if keyword == 'parameter':
free_param = _free_parameter_from_conf_record(name, value)
elif _VAR_DECL.search(keyword):
free_param = _free_parameter_from_var_line(name, keyword, value)
else:
continue
free_params.append(_with_start_point(free_param, start_points))
if not free_params:
raise PybnfError(
"No exportable free parameters found in the config (expected a 'parameter:' "
f"record or one of {sorted(EXPORTABLE_PRIOR_KEYWORDS)}).")
if start_points:
# A start point for a name no exportable free parameter claims. Silently dropping
# it is the failure this whole path exists to remove, and config.py refuses the
# same thing when the job is run.
names = ', '.join(sorted(start_points))
known = ', '.join(sorted(fp.name for fp in free_params))
raise PybnfError(
f"start point for unknown parameter(s) {names}",
f"The config declares a start point for {names}, which no exportable free "
f"parameter declaration names. The free parameters being exported are: "
f"{known}.")
return free_params
def _free_parameter_from_var_line(name, keyword, value):
"""The legacy positional ``<family>_var = <id> p1 [p2] [b|u]`` declaration."""
_require_exportable_prior(name, keyword)
# p1/p2 are the family's governing values (bounds for the Uniform families,
# loc/scale or shape/scale for the two-parameter location families). A one-parameter
# unbounded family (exponential/chisquare/rayleigh, #417) carries only p1.
p2 = float(value[1]) if len(value) >= 2 else None
# The 3rd token is the native reflecting-bounds flag, which only the Uniform families
# take ('u' -> False, 'b'/absent -> True). It used to be read as inert and dropped, so
# the exporter built a BOUNDED parameter from a conf that declared an unbounded search
# and wrote the box into PEtab's hard bounds without a word (#736). Passed through, so
# the parameter matches the one the fitter builds and petab_parameter_row can refuse
# the shape PEtab has no field for.
bounded = value[2] if len(value) >= 3 else True
return FreeParameter(name, keyword, float(value[0]), p2, bounded=bounded)
def _free_parameter_from_conf_record(name, fields):
"""A new-era ``parameter:`` record (ADR-0043) -> its ``FreeParameter`` (#733).
Built through :func:`~pybnf.parameter_record.free_parameter_from_record`, the same
mapping the fitter loads a job with, so the exported row describes the parameter the
fit would actually search rather than a second reading of the grammar. The record
resolves to one of the ordinary ``*_var`` keywords, so the exportability gate below is
the one the positional line goes through -- a record's boundaries are the same
boundaries, reached by a different spelling.
``initialization_distribution`` is fixed at ``prior`` rather than read from the config:
it selects where an algorithm draws its start points, which is run recipe rather than
problem, has no home in a PEtab table, and is what the positional line above defaults
to.
"""
try:
free_param = free_parameter_from_record(name, fields, INITIALIZATION_PRIOR)
except OutOfBoundsException:
# An out-of-box initial_value, as at the Configuration loader's own call site: the
# bare OutOfBoundsException subclasses Exception, so pybnf.main reports it as "an
# unknown error ... please report this bug" on a config the user wrote (#583).
raise PybnfError(
f"start point out of bounds for '{name}'",
f"Parameter '{name}' declares an initial_value outside its own lower/upper "
f"bounds. A start point is refused rather than moved into the box.",
"Correct the initial_value, or widen the parameter's bounds.")
_require_exportable_prior(name, free_param.type)
return free_param
def _require_exportable_prior(name, keyword):
"""Refuse a free parameter whose prior family PEtab v2 cannot state (ADR-0025, #423).
Keyed on the ``*_var`` keyword, which both declaration spellings resolve to, so a
``parameter:`` record hits the identical boundary as the positional line that builds
the same parameter.
"""
if keyword in EXPORTABLE_PRIOR_KEYWORDS:
return
if keyword in _NO_PRIOR_KEYWORDS:
raise NotImplementedError(
f"Free parameter '{name}' is a no-prior point start (a '{keyword}' line, or a "
f"'parameter:' record with an 'initial_value:' but no 'prior:' and no "
f"'lower:'/'upper:' box). A flat improper prior is not a PEtab probability "
f"family, and a PEtab estimated parameter needs bounds or a prior. Give it a "
f"prior or a box; the exporter writes {sorted(EXPORTABLE_PRIOR_KEYWORDS)} "
f"(ADR-0025, #423).")
raise NotImplementedError(
f"Free parameter '{name}' is a '{keyword}'; the exporter writes the PEtab prior "
f"families {sorted(EXPORTABLE_PRIOR_KEYWORDS)}, and PEtab v2 has no "
f"priorDistribution spelling for this one -- it defines no log- form for "
f"cauchy/gamma/exponential/chisquare/rayleigh, no natural-log ('ln') sampling "
f"scale, and no three-parameter family such as student_t (ADR-0025, #423).")
def _with_start_point(free_param, start_points):
"""Attach this parameter's ``start_point =`` line, if the config declares one (#719).
The merge rule is ``Configuration._load_start_point``'s (ADR-0117): a ``parameter:``
record's ``initial_value:`` and a ``start_point`` line are two spellings of one fact,
so they may both be present when they agree and are refused when they disagree --
silently preferring one would reintroduce the class of failure the start-point work
exists to remove.
"""
if free_param.name not in start_points:
return free_param
start = start_points.pop(free_param.name)
if free_param.value is not None:
if free_param.value != start:
raise PybnfError(
f"contradictory start point for '{free_param.name}'",
f"Parameter '{free_param.name}' is given two different start points: "
f"initial_value: {free_param.value} on its 'parameter:' record, and "
f"'start_point = {free_param.name} {start}'. Delete one of them.")
return free_param
try:
return free_param.set_value(start, reflect=False)
except OutOfBoundsException:
# A start point outside the parameter's own box. PEtab cannot express one either
# (a nominalValue outside lowerBound/upperBound is what the importer refuses on
# the way back in), and the bare OutOfBoundsException would reach the user as "an
# unknown error ... please report this bug" (#583).
raise PybnfError(
f"start point out of bounds for '{free_param.name}'",
f"The config starts '{free_param.name}' at {start}, which is outside the box "
f"[{free_param.lower_bound}, {free_param.upper_bound}] its own declaration "
f"gives it. A start point is refused rather than moved.",
"Correct the start point, or widen the parameter's bounds.")
def _start_points_from_conf(conf):
"""``{name: theta}`` for every ``start_point = <parameter> <value>`` line (#583).
One of the two spellings ``Configuration._load_start_point`` merges; the other, a
``parameter:`` record's ``initial_value:`` field, arrives on the built
``FreeParameter.value`` and is reconciled with this one in :func:`_with_start_point`.
"""
return {key[1]: float(value) for key, value in conf.items()
if isinstance(key, tuple) and len(key) == 2 and key[0] == 'start_point'}
@dataclass(frozen=True)
class _SbmlModelView:
"""An SBML model presented through the attribute surface the exporter's shared
classification reads (the :class:`~pybnf.petab._bngl.BnglEntities` subset), with SBML
semantics (ADR-0040 -- the export mirror of the ADR-0036 import dispatch):
* ``observable_names`` are the SBML **species** (the trajectory's bare-name output
columns). SBML has no BNGL-style observables, so a bare-name ``observableFormula``
names a species directly -- the inverse of the importer mapping a bare-name formula
back to a species column.
* ``function_names`` / ``function_bodies`` are **empty**: SBML has no global BNGL
functions, so an SBML observable is never inlined; an expression observable is carried
in the conf measurement-model layer (``observable: <id>, formula: <expr>``) and
classified via ``measurement_models``.
* ``parameters`` maps each global parameter id to its numeric nominal (or ``None``) --
the free-parameter binding set, the condition-target set, and the nominal source.
* ``compartment_names`` are the SBML compartments (also a valid condition target).
"""
text: str
parameters: dict
observable_names: frozenset
compartment_names: frozenset
function_names: frozenset = frozenset()
function_bodies: dict = field(default_factory=dict)
def _model_language(model_file):
"""The PEtab model language for a model file, by extension: ``.bngl`` -> ``'bngl'``,
``.xml`` -> ``'sbml'``. Mirrors ``config.py``'s model-suffix dispatch (which also maps
``.xml`` to its SBML backends)."""
if model_file.endswith('.bngl'):
return 'bngl'
if model_file.endswith('.xml'):
return 'sbml'
raise NotImplementedError(
f"Model '{model_file}' has an unrecognized extension; the exporter emits BNGL "
f"('.bngl') and SBML ('.xml') models (ADR-0025/0040). An Antimony ('.ant') model "
f"would need an Antimony->SBML conversion first -- a later chunk.")
def _read_model(model_file, path, language):
"""Read the job's model into the entity view the exporter consumes, per ``language``.
BNGL -> :class:`~pybnf.petab._bngl.BnglEntities`; SBML -> :class:`_SbmlModelView`
(the same attribute surface, SBML semantics). Both expose ``text`` (verbatim source),
``parameters`` (the bindable ids + nominals), ``observable_names`` (bare-name observable
columns), ``function_names``/``function_bodies`` (BNGL only), and ``compartment_names``
-- everything the language-agnostic classification reads. The dispatch is the export
peer of the importer's :func:`~pybnf.petab.import_._model_namespace`.
"""
text = Path(path).read_text(encoding='utf-8', errors='replace')
if language == 'sbml':
ent = parse_sbml_model(text)
return _SbmlModelView(
text=text,
# Every global parameter id -> its nominal (float or None); the binding set is
# the key set, the nominal source is the value (a value-less parameter -> None).
parameters={p: ent.parameter_values.get(p) for p in ent.parameter_names},
observable_names=ent.species_names,
compartment_names=ent.compartment_names)
return parse_bngl_model(text)
# ---------------------------------------------------------------------------
# Emitting the PEtab-clean model and problem.yaml
# ---------------------------------------------------------------------------
# The hand-written actions the exporter may drop, by BNGL action name (#900). Under edition 2
# the fitter runs every action the model file holds -- in ``begin actions`` or loose after
# ``end model`` -- except ``generate_network`` and ``setOption`` (``BNGLModel.__init__``), and
# then the simulations it builds from the ``experiment:`` lines (``BNGLModel.add_action``),
# each network-based one preceded by ``resetConcentrations()``. PEtab has no such preamble, so
# an action may be dropped only if it cannot change what those simulations start from. Checked
# by reading ``pset.py`` and BioNetGen 2.9.3's ``Perl2`` sources, and by running the same fit
# with and without each action under BNG2.pl and under bngsim:
#
# * a simulation changes species amounts only, and the reset before each network-based
# experiment restores the seed (or the last ``saveConcentrations()`` snapshot, and saving is
# refused). A network-free (``method: nf``) experiment gets no reset, so it continues from
# where a hand-written simulation left off; that pairing is refused separately.
# ``method=>"protocol"`` runs the model's ``begin protocol`` block, whose own actions can set
# parameters, so it is not a plain simulation and is refused.
# * ``resetConcentrations()`` restores the seed or a saved snapshot, and saving is refused.
# * the ``write*`` / ``visualize`` actions write files and change nothing.
#
# Everything else is refused, the unknown included. ``parameter_scan`` / ``bifurcate`` look
# like simulations but are not droppable: BNG2.pl leaves the scanned parameter at its last scan
# value (``BNGAction.pm``'s scan never restores it), which every later experiment then sees
# (measured: the fit's objective moved from 7e-12 to 1606), while bngsim restores it -- so the
# fitter itself disagrees by backend, and no export can match both.
_SIMULATION_ACTIONS = frozenset({
'simulate', 'simulate_ode', 'simulate_ssa', 'simulate_pla', 'simulate_psa', 'simulate_nf'})
_DROPPABLE_ACTIONS = frozenset({
'resetConcentrations',
'writeXML', 'writeSBML', 'writeNetwork', 'writeNET', 'writeFile', 'writeModel', 'writeBNGL',
'writeMfile', 'writeMexfile', 'writeMEXfile', 'writeMDL', 'writeLatex', 'writeSSC',
'writeSSCcfg', 'writeCPPfile', 'writeCPYfile', 'visualize'})
# The network-free experiment methods (``Action.VALID_METHODS``). ``BNGLModel.add_action``
# omits the reset for ``nf`` alone today; RuleMonkey (``rm``) runs on the same bngsim
# network-free session, which keeps its state from one simulate to the next and does not run a
# reset, so it is counted too rather than trusting a reset line that session skips.
_NETWORK_FREE_METHODS = frozenset({'nf', 'rm', 'rulemonkey'})
# An action's name: the identifier before its opening parenthesis.
_ACTION_NAME = re.compile(r'([A-Za-z_]\w*)\s*\(')
_PROTOCOL_METHOD = re.compile(r'method\s*=>\s*["\']protocol["\']')
# The line ``BNGLModel._synthesized_generate_network_line`` writes when no ``generate_network``
# conf key is set: the default network generation every BNGL consumer applies unasked, so it
# is left implicit in the exported model rather than written into every one of them.
_BARE_GENERATE_NETWORK = 'generate_network({overwrite=>1})'
def _action_code(raw):
"""One ``BNGLModel.actions`` entry as code: comments stripped, the physical lines of a
backslash continuation joined, whitespace collapsed. Empty for a blank or comment line."""
code = re.sub(r'#[^\n]*', '', raw)
code = re.sub(r'\\\s*\n', ' ', code)
return ' '.join(code.split())
def _require_droppable_actions(model, model_file, experiment_methods):
"""Refuse a BNGL model whose hand-written actions the export cannot drop without changing
the fit (#900); see ``_DROPPABLE_ACTIONS`` for which actions are droppable and why.
``model.actions`` is the fitter's own list -- the scan ``BNGLModel.__init__`` does for the
fit -- so this checks exactly the lines the fit runs ahead of its experiments, the ones in
``begin actions`` and the loose ones after ``end model`` alike. ``experiment_methods`` maps
the model's experiments to their ``method:`` (``None``: unknown, none taken as
network-free)."""
refused, simulations = [], []
for raw in model.actions:
code = _action_code(raw)
if not code:
continue # a blank or comment line inside the block
match = _ACTION_NAME.match(code)
name = match.group(1) if match else None
if name in _SIMULATION_ACTIONS and not _PROTOCOL_METHOD.search(code):
simulations.append(code)
elif name not in _DROPPABLE_ACTIONS:
refused.append(code)
if refused:
raise NotImplementedError(
f"Model '{model_file}' carries action(s) that the fit runs before its experiments "
f"but the PEtab export would drop: {'; '.join(refused)}. Under edition 2 the fitter "
f"runs every action in the model file (in 'begin actions' or loose after 'end "
f"model'), except generate_network and setOption, ahead of the simulations it builds "
f"from the 'experiment:' lines, so such an action can change what every experiment "
f"starts from -- a parameter value, a species amount, or the snapshot "
f"resetConcentrations() restores (a parameter_scan or bifurcate leaves its parameter "
f"at the last scanned value under BioNetGen). PEtab has no step to carry it, so the "
f"exported problem would have a different optimum. Move the change into the model's "
f"parameters or seed species, or into a 'condition:' on the experiments, or delete "
f"the line if the fit should not use it. The export drops only simulate*, "
f"resetConcentrations and write*/visualize, and keeps generate_network (#900).")
network_free = sorted(name for name, method in (experiment_methods or {}).items()
if method in _NETWORK_FREE_METHODS)
if simulations and network_free:
raise NotImplementedError(
f"Model '{model_file}' carries hand-written simulation action(s) "
f"({'; '.join(simulations)}) that the fit runs ahead of experiment(s) "
f"{network_free}, which run network-free (method: nf / rm). The fitter starts a "
f"network-free experiment without resetting the species, so it continues from the "
f"state the hand-written simulation left, while PEtab starts every experiment from "
f"the model's initial state. Delete the hand-written simulation action(s) from the "
f"model to export (#900).")
def _exported_generate_network_line(model, text_lines, experiment_methods):
"""The ``generate_network`` line the exported model carries: the fitter's own (#485/#901).
* The model has its own line: the fitter uses it and ignores the ``generate_network`` conf
key ("an explicit line in the model always wins", ``pset.py``) -- the last one if there
are several, since the scan keeps overwriting ``generate_network_line``. So does the
export.
* Otherwise the fitter synthesizes ``generate_network({overwrite=>1,<opts>})`` from the key
whenever it generates a network -- for a hand-written network-based simulation, or for
any experiment that is not network-free (``BNGLModel.add_action``). The export used to
write nothing, so a cap stated in the job was lost and a PEtab consumer built a
different, or an unbounded, network (#901). The synthesized line is now written; the bare
default (no key) is what every consumer does unasked, so it stays implicit.
``None`` when the exported model needs no line."""
own = any(re.match('generate_network', text_lines[i].split('#', 1)[0].strip())
for i in model.action_line_indices)
if own:
return model.generate_network_line.strip()
generates = (model.generates_network or experiment_methods is None
or any(method != 'nf' for method in experiment_methods.values()))
if not generates:
return None
line = model._synthesized_generate_network_line()
return None if line == _BARE_GENERATE_NETWORK else line
[docs]
def clean_model_for_petab(text, model_file='model.bngl', generate_network_options=None,
experiment_methods=None):
"""Return a PEtab-clean copy of a BNGL model: its actions removed, and the network
definition the fit used (``generate_network``) kept.
New-era BNGL binds free parameters **by id** (ADR-0034), so the source model already
carries bare parameter ids with real nominal values -- exactly what PEtab estimates.
"PEtab-clean" therefore drops the *simulation* actions (PEtab drives simulation via the
measurement times / experiments, not the model's own ``simulate`` calls) while keeping
``generate_network`` -- a network-definition / compilation directive, not a simulation
action. That directive carries the model's finiteness cap (``max_stoich`` / ``max_agg`` /
``max_iter``); dropping it would silently turn a model that is finite only under the cap
into one that network-generates unbounded, with no error or warning (#485). The line
written is the one the fitter runs (:func:`_exported_generate_network_line`): the model's
own, or the one the fitter synthesizes from the job's ``generate_network`` key (#901). The
exported model then carries its own cap, so ``import_.py`` (which copies the model
byte-verbatim) round-trips it and any BNG2.pl / PyBNF consumer stays finite; with no line
to keep, the actions disappear entirely.
The model is read with the fitter's own scanner (``BNGLModel``), so the lines removed are
exactly the ones the fit reads as actions -- the ``begin actions`` block and any loose
action after ``end model`` -- and the kept line goes in a minimal ``begin actions`` block
at the end of the file, where the fitter writes it. ``setOption`` and its siblings stay
where they are (the fitter keeps them in the model text too), as do a comment outside the
actions block and a protocol block. An action the export cannot drop without changing the
fit -- one that sets a parameter or species, saves a snapshot, or reads a file -- raises
rather than vanishing (:func:`_require_droppable_actions`, #900). The reaction network and
the ``begin functions`` block -- which carry the measurement model -- are carried verbatim.
A fit-and-mutated parameter keeps its model name (``v1``) here as a plain nominal-valued
parameter (always overridden by its Condition); only the parameter *table* carries the
surrogate ``v1__REF`` (ADR-0027).
``generate_network_options`` is the job's ``generate_network`` key and
``experiment_methods`` maps the experiments on this model to their ``method:``
(:func:`_experiment_methods_by_model`); the ``None`` defaults mean no key and network-based
experiments. ``model_file`` names the model in error messages.
A legacy ``<name>__FREE`` marker in the model **code** is **rejected**: new-era binds by
id, so a model still carrying one was not modernized, and shipping it would dangle an
undefined ``v1__FREE`` symbol in PEtab. The error names the bind-by-id contract rather
than letting the PEtab oracle reject it opaquely. The scan ignores ``#`` line comments --
a comment may legitimately mention the retired ``KD1__FREE`` form as a counter-example
(``receptor_v2.bngl`` does), which is documentation, not a dangling binding.
"""
code = re.sub(r'#[^\n]*', '', text) # strip line comments before the marker scan
if _FREE_TOKEN.search(code):
raise PybnfError(
"This BNGL model carries a legacy '__FREE' marker, but PEtab export is a "
"new-era feature where free parameters bind by id (ADR-0034). Declare the "
"model's fit parameters as bare ids with nominal values (e.g. 'v1 0.5', not "
"'v1 v1__FREE') and list them as free parameters in the .conf.")
try:
model = BNGLModel(model_file, suppress_free_param_error=True,
generate_network_options=generate_network_options, text=text)
except ModelError as exc:
raise PybnfError(f"Model '{model_file}' could not be read as BNGL: {exc}.") from exc
_require_droppable_actions(model, model_file, experiment_methods)
lines = text.splitlines(keepends=True)
network_line = _exported_generate_network_line(model, lines, experiment_methods)
out = ''.join(line for i, line in enumerate(lines) if i not in model.action_line_indices)
if network_line is not None:
if out and not out.endswith('\n'):
out += '\n'
out += f'begin actions\n{network_line}\nend actions\n'
return out
[docs]
def write_problem_yaml(path, models, has_conditions=False, has_experiments=False,
has_mapping=False):
"""Write a PEtab v2 ``problem.yaml`` referencing the tables and the model(s).
``models`` is a list of ``(model_id, location, language)`` tuples (ADR-0041): one entry
for a single-model job (byte-identical to the pre-multi-model output), or N entries in
declaration order for a multi-model job. ``language`` is each model's PEtab language
(``bngl`` or ``sbml``, ADR-0040) -- emitted verbatim so every model declares its own
native language (a BNGL + SBML mix is just two entries). ``has_mapping`` adds the
``mapping_files`` entry for the species-amount mapping table (ADR-0062)."""
parts = [
'format_version: 2.0.0\n',
'parameter_files:\n - parameters.tsv\n',
'observable_files:\n - observables.tsv\n',
'measurement_files:\n - measurements.tsv\n',
]
if has_conditions:
parts.append('condition_files:\n - conditions.tsv\n')
if has_experiments:
parts.append('experiment_files:\n - experiments.tsv\n')
if has_mapping:
parts.append('mapping_files:\n - mapping.tsv\n')
parts.append('model_files:\n')
for model_id, location, language in models:
parts += [
f' {model_id}:\n',
f' location: {location}\n',
f' language: {language}\n',
]
Path(path).write_text(''.join(parts))
