PyBNF gradient plumbing (pybnf.gradient)

The pybnf.gradient package holds the PyBNF-side machinery that turns a simulation backend’s forward output-sensitivity tensor into objective gradients and residual Jacobians for the gradient-based optimizers (trf/lbfgs). It has three parts: a routing layer that maps edition-2 free parameters onto each experiment’s sensitivity request and per-parameter chain-rule factor (pure mapping, no objective math); an assembly layer that combines the sensitivity tensor with that routing into an objective gradient / residual Jacobian in sampling space (the form the optimizer consumes); and an errors module whose GradientNotSupported signals a configuration outside the differentiable set, so a caller can fall back to a gradient-free step.

The capability gate and the per-layer math (fixed and estimated noise scale, log/lognormal scale, per-observable transforms and normalization, the Laplace and Student-t families with mean centering, and constraint-penalty gradients) are documented for users in Gradient-based fitting (forward sensitivities).

Errors

Gradient-path exceptions (#449/#385).

A single, dependency-free home for the capability-gate exception so both the objective seam (LikelihoodObjective.residual_point) and the assembly (pybnf.gradient.assembly) can raise/catch it without an import cycle: objective.py -> gradient.errors is one-way (the assembly takes the objective as an argument and never imports it back).

exception pybnf.gradient.errors.GradientNotSupported(log_message, user_message=None)

The gradient path cannot (yet) differentiate this objective configuration.

Cut-1 of the #385 epic assembles the gradient for exactly one case: the default Gaussian family, additive on the LINEAR scale, with the prediction as the MEDIAN and a fixed (non-estimated) sigma. Every other configuration – estimated sigma (layer D), a log/lognormal scale (layer E), a trajectory transform (layer F), an asymmetric family (layer G), a constraint penalty (layer I), … – raises this, naming what is unsupported, until its layer lands. It is a PybnfError so a caller (#386’s optimizer) can fall back to a gradient-free step with a clean, actionable message instead of a silently wrong derivative.

Routing

Free-parameter -> forward-sensitivity routing for the gradient path (#448, #385).

Step B of the #385 gradient-plumbing epic: a pure mapping, no objective math. Given an edition-2 fit’s free parameters, a model’s bind-by-id namespace (ADR-0034), and an experiment’s condition perturbation (ADR-0028), this computes – per experiment – which free parameters become a bngsim forward-sensitivity request and, for each, the chain-rule factor that converts a native sensitivity column into the derivative w.r.t. the free parameter.

The routing feeds two consumers:

• the sensitivity_params / sensitivity_ic lists handed to #447’s request at each experiment’s Simulator (pybnf.bngsim_model.net_model.BngsimModel.enable_output_sensitivities()); and

• the per-free-parameter factor that #449 multiplies into the per-experiment objective Jacobian. Building that Jacobian is out of scope here.

Bind by id (ADR-0034)

A parameter: free parameter (ADR-0043) binds to the model entity of the same id – the contract set_param already uses. So classification is by name (classify_free_param()):

• a free parameter whose id is a species’ bare initial-value expression (_parse_net_species_initializers) routes to sensitivity_ic keyed by the species. An initial-condition parameter does not appear in the ODE RHS, so the parameter sensitivity axis would be identically zero; bngsim’s ic axis carries it. (Checked first, because an IC parameter is also a begin parameters id.)

• a free parameter whose id is a model parameter routes to sensitivity_params keyed by the parameter.

• a free parameter matching no model id (e.g. a free sigma) is not a sensitivity request – it carries no model column (its gradient is assembled in layer D, #449+). This is the existing unmatched-parameter path (a warning, not an error – ADR-0034).

Per-condition perturbation = local derivative (ADR-0028)

A condition: is a pybnf.pset.MutationSet on the base model; each pybnf.pset.Mutation carries operation in {=,+,-,*,/}. For a free parameter p perturbed in an experiment’s condition the chain-rule factor is the local derivative d(perturbed p)/dp (condition_factor()):

• p = c -> p is pinned to a constant in that experiment => factor 0 (its column is dropped from the request);

• p * c -> factor c; p / c -> factor 1/c;

• p + c / p - c -> factor 1 (a shift has unit slope);

• an unperturbed free parameter -> factor 1.

Multiple mutations on the same id compose: affine maps compose, so the multiplicative parts multiply; any = pins the parameter, driving the composed factor to 0.

Cut-1 scope

IC routing matches a bare initializer (species <- p). A parameter that both appears in the ODE RHS and seeds an initial value, or a species seeded by a non-trivial expression of the free parameter (2*p), is a later layer – the clean param/ic split is what #449’s first FD-oracle case needs.

class pybnf.gradient.routing.ExperimentRouting(routes: dict)

The per-experiment routing object: {free_param -> ParamRoute} plus the derived sensitivity_params / sensitivity_ic request lists handed to #447’s gradient path.

property sensitivity_ic

sensitivity_ic= for the experiment’s Simulator: the species of every IC-bound free parameter with a non-zero factor, de-duplicated in declared free-parameter order.

property sensitivity_params

sensitivity_params= for the experiment’s Simulator: every parameter-bound free parameter with a non-zero factor (a pinned = column is dropped), de-duplicated in declared free-parameter order.

class pybnf.gradient.routing.ParamRoute(free_param: str, target: str, key: object, factor: float)

How one free parameter maps onto an experiment’s forward-sensitivity request.

target is PARAM (kinetic/global -> sensitivity_params), IC (species initial value -> sensitivity_ic), or NONE (bound to no model id – a free sigma; no model column). key is the request key the tensor is read by: the parameter id for PARAM, the species for IC, None for NONE. factor is the chain-rule derivative d(perturbed param)/d(free param) for this experiment’s condition (#449 multiplies it into the Jacobian column); a pinned (=) parameter has factor 0 and is dropped from the request list.

pybnf.gradient.routing.apply_routing(model, routing)

Hand a routing’s request lists to #447’s gradient path on model.

Calls BngsimModel.enable_output_sensitivities() with the routing’s sensitivity_params / sensitivity_ic – the capability-gated activation of the gradient path. A build without forward output sensitivities refuses there (#447). Returns the same routing for chaining.

pybnf.gradient.routing.classify_free_param(free_param, param_ids, species_initializers)

Classify one free parameter by id (ADR-0034): return (target, key).

Checks the species initial-value namespace first: a free parameter that is a species’ bare initial-value expression routes to the IC axis keyed by the species (an IC parameter is absent from the ODE RHS, so the parameter axis is identically zero). Otherwise a match in the begin parameters namespace routes to PARAM keyed by the id; no match at all is NONE (a nuisance such as a free sigma – no model column).

param_ids is the model’s begin parameters namespace (any container supporting in); species_initializers is the (species, initial-expr) list from _parse_net_species_initializers.

pybnf.gradient.routing.condition_factor(free_param, condition)

The chain-rule factor d(perturbed param)/d(free param) for one free parameter under an experiment’s condition (ADR-0028).

condition is the experiment’s pybnf.pset.MutationSet, or None for the unperturbed wildtype. Composes every mutation that targets this id (affine maps compose, so the multiplicative parts multiply): = contributes 0 (pins the parameter to a constant), *c contributes c, /c contributes 1/c, and + / - contribute 1 (an additive shift has unit slope). An id the condition does not touch keeps the identity factor 1.

pybnf.gradient.routing.route_experiment(free_params, param_ids, species_initializers, condition=None)

Build the ExperimentRouting for one experiment (pure – no model, no sim).

free_params is the ordered free-parameter id list (the config’s declared variables); param_ids the model’s begin parameters namespace; species_initializers the (species, initial-expr) pairs; condition the experiment’s pybnf.pset.MutationSet (None for the wildtype experiment).

pybnf.gradient.routing.route_for_model(model, free_params, condition=None)

route_experiment() against a live model’s bind-by-id namespaces.

Reads the model’s begin parameters ids and (species, initial-expr) pairs through BngsimModel.sensitivity_entity_namespace() (the only model coupling), so the routing core stays backend-agnostic. condition may be a pybnf.pset.MutationSet, a condition name resolved against model.mutants, or None for the wildtype experiment.

Assembly

Gaussian objective gradient + residual-Jacobian assembly (#449, #385).

Step C of the #385 gradient-plumbing epic. Given #447’s per-experiment forward output-sensitivity tensor (Data.output_sensitivities) and #448’s per-experiment routing (pybnf.gradient.routing.ExperimentRouting), assemble – for the default Gaussian, LINEAR-scale, fixed-sigma objective, summed across experiments – both forms of the gradient:

1. the scalar dF/du (for quasi-Newton / L-BFGS-B), and

2. the residual vector + residual-Jacobian (for trust-region least-squares, #386’s primary path).

Convention pin (issue #449)

PyBNF’s Gaussian loss is data_fit = (pred - obs)**2/(2 sigma**2) with mu = pred for the default (LINEAR, MEDIAN) family, so per scored point i:

• residual rho_i = (pred_i - obs_i)/sigma_i, loss = 1/2 ||rho||**2;

• residual-Jacobian (native param space) J_ij = (1/sigma_i) * factor_j * d pred_i/d theta_j – d pred/d theta from #447’s tensor, factor_j from #448’s routing;

• scalar gradient dF/d theta = J^T rho (not 2 J^T rho).

scipy.least_squares minimizes 1/2 ||rho||**2 with the same rho/J, so the residual form and the scalar form agree by construction – the optimizer walks the surface PyBNF reports. The per-point bootstrap weight w_i (1.0 unless bootstrapping) is folded in as sqrt(w_i) on both rho_i and J_i, so 1/2 ||rho||**2 stays sum_i w_i * data_fit_i – exactly what evaluate sums (eval_point * weight).

Native -> sampling space (once, ADR-0029)

rho is scale-invariant; the Jacobian moves to the sampling space the optimizer walks by J -> J @ diag(d theta/d u), applied once at the end. d theta/d u is one autodiff of each parameter’s scale inverse_jax (priors/scale.py) – no hand-written per-scale derivative. A LINEAR parameter has d theta/d u = 1 and is short-circuited, so the common (all-linear) path needs no jax; only a log-scaled parameter pulls in the optional pybnf[jax] extra (the house pattern, ADR-0019).

Estimated noise scale (layer D, #451)

An estimated sigma – the edition-2 noise_model = normal, sigma = fit <param> surface (ADR-0021/0034), a freely-named free parameter; equivalently chi_sq_dynamic’s legacy sigma__FREE default – keeps the Gaussian normalizer, so the per-point loss is (pred-obs)**2/(2 sigma**2) + log sigma and gains a sigma column d loss/d sigma = -(pred-obs)**2/sigma**3 + 1/sigma (objective.noise_grad_point). The routing binds that free parameter by id (ADR-0034); estimated noise is matched by source type (FreeParameterSigma), never by a name convention. +log sigma is not a sum of squares, so it cannot live in the residual/Jacobian form: this assembly adds the sigma column straight to the scalar gradient and leaves the residual-Jacobian a faithful least-squares model of the data fit alone – flagged by GradientResult.least_squares_exact (False once any estimated scale is present), so #386’s trust-region path knows to use the scalar gradient (L-BFGS) for an estimated-sigma fit. The free sigma routes to NONE in #448 (no model column), so its gradient comes entirely from this normalizer + the sigma-dependence of the data fit, never from the sensitivity tensor.

Trajectory transforms + normalization (layer F, #453)

_prediction may form the scored value from the raw observable through a per-observable transform: a cumulative -> incident difference (ADR-0051), a per-measurement scale/ offset formula (ADR-0045), or an upstream Data-level normalization (ADR-0053). Each makes ∂pred/∂θ differ from the raw observable sensitivity, so the assembly reads a raw_sens(col, row) accessor (the #447 tensor, routing-factor-folded, with normalization’s own quotient/chain rule threaded in – _normalized_sensitivity) and hands it to the objective’s prediction_sensitivity() seam, which mirrors _prediction branch for branch (cumulative differences sensitivity rows; a per- measurement formula chains its symbolic gradient through each referenced column’s sensitivity plus any estimated placeholder it names – which, unlike a free σ, does enter ∂pred/∂θ and so lands in the residual-Jacobian). A plain column collapses to the raw sensitivity, so the no-transform path is byte-identical.

Asymmetric / non-Gaussian families (layer G, #454; least-squares residual #459)

A family yields an exact least-squares residual/Jacobian iff its data fit reformulates as a smooth half-square. Two do: the Gaussian (data_fit = 1/2 rho**2) and – the #459 follow-up – the Student-t, whose exact square-root-loss residual r = sign(z) sqrt(2 data_fit) = sign(z) sqrt((ν+1) log1p(z²/ν)) satisfies both 1/2 r² == data_fit and r · d r/d pred == d(data_fit)/d pred and is smooth through z=0 (r ~ sqrt((ν+1)/ν) z, downweighting the tails). Both route through residual_point (NoiseModel.residual / d_residual_d_prediction), contributing a residual row + native Jacobian row, so a fixed-scale Student-t fit is least_squares_exact – the Gaussian’s exact-least-squares status recovered for the robust family, and #386’s LM/TRF can fit it directly.

Laplace has no such clean residual: its L1 data fit |·|/b gives sqrt(2 data_fit) ~ sqrt|z|, a cusp with infinite slope at z=0, so it stays scalar-only (the count family likewise). Such a family’s objective gradient is assembled from the universal scalar form sum_i w_i · d(data_fit_i)/d(pred_i) · d(pred_i)/d θ – the per-family slope d(data_fit)/d(pred) (objective.data_fit_grad_point -> NoiseModel.d_data_fit_d_prediction) chained through the same layer-F prediction_sensitivity d pred/d θ, accumulated into a separate data_fit_gradient vector exactly as the estimated-noise column accumulates into noise_gradient; the result is flagged not least_squares_exact (the residual/Jacobian then model only the Gaussian/Student-t columns, if any), the signal #386’s trust-region path uses the scalar gradient (L-BFGS). The estimated-noise column generalizes per family throughout: Laplace’s b, and Student-t’s sigma and df (the first multi-parameter estimated-noise gradient, ADR-0058) – a retained normalizer is never a square, so an estimated-scale fit is inexact for any family (its data-fit residual, if it has one, still stacks; the normalizer column rides noise_gradient). The prediction (the MEDIAN or, layer G, the MEAN) and the noise scale (linear or log, layer E) compose throughout; objective.has_least_squares_residual routes each column. An all-Gaussian fit never touches data_fit_gradient (it stays zero), so that path is byte-identical. The negative-binomial count family rides the same scalar data_fit_gradient path (#458): its prediction slope is the NB score chained through the median CDF-inversion implicit derivative, and its estimated dispersion r – self-normalizing PMF, so the whole column is in the data fit – the noise_gradient path. The dispersion column closes for MEAN and MEDIAN alike: a MEDIAN mean is solved from r, so its column folds in d mean/d r (the betainc first-parameter implicit derivative), the count analogue of the mean-on-log offset coupling the location-scale families fold in (#385/#458).

Constraint / qualitative penalties (layer I, #456)

A fit may add qualitative / inequality constraints (BPSL .con / .prop files) whose penalty is added to the objective. assemble_constraint_gradient() is the sibling assembler: each constraint’s penalty is a piecewise (static) or Gaussian-CDF (likelihood) function of an at-/ between-time readout q1 - q2, so its gradient is that readout’s forward sensitivity (read via a (model, suffix, observable)-keyed accessor over the #447 tensor + #448 routing) times the local penalty slope (penalty_gradient()). Like an estimated- noise normalizer, a penalty is not a sum of squares, so it lives on the scalar gradient only: a fit with active constraints is not least_squares_exact, and #386 adds this term to the objective gradient.

Measurement-model layer (layer H, #455)

A scored observable may be a materialized measurement-model column (ADR-0036): an expression observableFormula the objective’s MeasurementLayer adds to the trajectory before scoring (the SBML/Antimony / new-era PEtab path, where the SBML backend exposes the same species: forward sensitivities the net backend does for observables). Such a column is not in the #447 tensor, so the raw_sens accessor recognizes it and delegates to the model’s prediction_sensitivity() – the formula’s exact chain rule over each referenced column’s sensitivity (read back through this same accessor, so routing / normalization fold in) plus any fit parameter the formula names directly (an observation-model scale/offset, which – like a per-measurement scale – enters ∂pred/∂θ and lands in the residual-Jacobian). A plain (non-measurement) column collapses to the tensor/normalized sensitivity, so the no-measurement path is byte-identical.

class pybnf.gradient.assembly.GradientResult(residual: ndarray, jacobian: ndarray, gradient: ndarray, param_names: list, least_squares_exact: bool = True)

The assembled gradient of a Gaussian objective at one parameter point.

residual is the stacked standardized residual rho (one entry per scored observation across all experiments, sqrt(weight)-folded). jacobian is the matching (n_obs, n_param) residual-Jacobian in sampling space (the d theta/d u transform already applied). gradient is the scalar dF/d u over the free parameters, in param_names order.

With a fixed-scale fit whose families all carry an exact least-squares residual – the Gaussian (any scale/location) and, #459, the Student-t (its smooth square-root-loss residual) – the data fit IS the whole objective, so the residual and scalar forms agree by construction (gradient == jacobian.T @ residual, 0.5||rho||**2 == evaluate) and least_squares_exact is True. It is False once the residual-Jacobian is no longer the whole story:

• an estimated noise parameter (layer D/G, #451/#454) – its retained normalizer (+log sigma, log(2 b), the df-block) is not a square, so it is folded into the scalar gradient only and the residual-Jacobian models the data fit alone; or

• a family with no clean least-squares residual (Laplace, whose L1 data fit is the cusp sqrt|z|; the count family, layer G #454/#459) – its data fit is not a smooth sum of squares, so it carries no residual at all and its whole data-fit gradient is on the scalar path.

Either way the scalar gradient is complete (jacobian.T @ residual over the residual-bearing columns – Gaussian / Student-t – if any, plus the data-fit and noise columns), and least_squares_exact = False is the signal that #386’s trust-region step must consume it rather than the bare residual.

pybnf.gradient.assembly.assemble_constraint_gradient(constraint_sets, sim_data_dict, routings, free_params)

The scalar gradient of the total constraint penalty w.r.t. the free parameters (layer I, #456/#385), in sampling space – the term #386 adds to the objective gradient for a fit with active constraints.

constraint_sets is an iterable of ConstraintSet; each constraint’s penalty is a piecewise (static) or Gaussian-CDF (likelihood) function of an at-/ between-time readout, so its gradient is that readout’s forward sensitivity times the local penalty slope (penalty_gradient()). sim_data_dict is the {model: {suffix: Data}} the penalties are scored on, each Data carrying the #447 sensitivity tensor; routings maps (model, suffix) -> ExperimentRouting (#448) so a readout’s sensitivity is factor-folded into the free-parameter axes exactly as the objective’s is. free_params fixes the parameter (column) order and supplies the native->sampling transform applied once at the end.

A penalty is not a sum of squares, so – like an estimated noise normalizer (layer D) – it lives on the scalar gradient only: a fit with active constraints is not least_squares_exact, and #386 must consume the scalar gradient. Returns a (n_param,) vector (zeros for an empty constraint set).

pybnf.gradient.assembly.assemble_gaussian_gradient(objective, experiments, free_params)

Assemble the scalar gradient and residual-Jacobian, summed across experiments.

objective is the fit’s LikelihoodObjective; it supplies each residual-bearing point’s residual through residual_point (a Gaussian or, #459, a Student-t – the families whose data fit is a smooth half-square), each no-residual (Laplace / count) point’s data-fit gradient through data_fit_grad_point (routed by has_least_squares_residual), and any estimated-noise gradient column through noise_grad_point (each gating the supported configuration – a Gaussian / Laplace / Student-t / negative-binomial family, MEDIAN or MEAN, any noise scale, noise fixed or single free parameters – raising GradientNotSupported otherwise). experiments is an iterable of (sim_data, exp_data, routing) triples – one per scored model/condition; each sim_data must carry the #447 output_sensitivities payload (the gradient path active), and routing is that experiment’s ExperimentRouting. free_params is the ordered list of FreeParameter defining the u-vector: it fixes the column order of the Jacobian and the entries of the scalar gradient, and supplies each parameter’s scale (current value -> d theta/d u).

Returns a GradientResult. The per-experiment routing is built once by the caller (#386) – it depends only on model structure, conditions, and free-parameter ids, never on the parameter values – so this per-evaluation assembly only reads the freshly simulated sensitivity tensors.