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, hint=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) may reach the trajectory on either axis, and which one is the backend’s answer, not an inference (ADR-0100, #537). d_param is the total derivative – the right-hand-side path plus whatever ∂x(0)/∂θ the backend seeded – so the rule is: route the id’s own sensitivity_params axis, and add a sensitivity_ic term only for a (species, param) pair BngsimModel.backend_ic_sensitivity() reports absent. Adding one the backend already carries doubles that column; omitting one it does not carry deletes the seeding. ode_rhs_symbols survives only as an optimization, dropping an axis that is provably identically zero (nothing seeded into it, and a right-hand side the model states never reads the id), so a model that cannot say simply keeps it.

• 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.

Reached through a condition (a per-condition estimated initial condition, ADR-0076)

A free parameter that binds no model id of its own can still reach the model through a condition that sets a model entity to its value – target = free_param (is_param_ref). The referenced free parameter then gains a RouteContribution on every sensitivity column target reaches (classify_condition_target()): PARAM on target itself for an ordinary global (e.g. a shared rate multiplier), plus one term per entity whose initial value target seeds – a species initial condition (IC) or a parameter an initialAssignment derives (PARAM). Because one free parameter may be assigned to several targets in one condition, and one target may seed several entities, a route is a sum over its contributions (ParamRoute).

The seed derivative is not assumed to be 1 (#530)

Each seeding term carries its own d(entity)/d(target), differentiated symbolically from the model’s initial-value expression by pybnf.gradient.derivative: 1 for the bare species = p, -1 for S_ = N_ - I0_, 2 for 2*p, and a compartment-unit factor where the species’ value and its assignment disagree on amount vs concentration. A derivative that is not a bare number – d(beta_N)/d(R0_) = gamma_/N_ – is point-dependent: it is carried symbolically and evaluated at each evaluated PSet by ExperimentRouting.at_point().

One contribution per native column (#537)

A route’s derivative is the sum over its contributions, so two contributions naming the same (target, key) column would have the assembly read one tensor column twice and add it twice. That is legitimate arithmetic when the two terms are genuinely different chain-rule paths that happen to meet – two condition targets seeding one species’ initial value – and it is indistinguishable, once assembled, from a column counted twice by mistake. So route_experiment() folds same-column terms into one contribution whose derivative tree is their sum, making “one contribution per native column” a structural property of every routing; ExperimentRouting.check_column_multiplicity() then lets the assembly assert that property per experiment rather than assume it (#537: a Raia_CancerResearch2011 gradient column came back at exactly twice its central difference, once, on the one route in that fit with an ic-axis contribution, and never again).

The fold on its own would defeat that check, though – it turns two same-column terms into one term of doubled factor, which after the fact is indistinguishable from a legitimate single term. So each contribution carries the chain-rule path(s) it came from (RouteContribution.origins: 'bind' for the bind-by-id classification, 'ref:<target>' per condition parameter-reference), and the check reads those. One path reaching one column twice is the defect; two paths meeting on it is arithmetic, and only the labels tell them apart once the factors have been summed.

Scope

The seed grammar is arithmetic (+ - * / **, numbers, symbols). An initial value reached through a function call, an assignmentRule, or a second initialAssignment is a chain this routing does not compose: classify_condition_target() raises GradientNotSupported rather than emit a wrong or silently-zero column. Bind-by-id routing of a free parameter that itself seeds several species remains single-column (classify_free_param()).

class pybnf.gradient.routing.ExperimentRouting(routes: dict, nominal_values: dict = <factory>, condition: object = None)

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

nominal_values (the model’s {parameter id: value} table) and condition are retained only so at_point() can rebuild the environment a point-dependent seed factor is evaluated in (#530); a routing whose factors are all constants ignores them.

at_point(param_values)

This routing with every point-dependent factor evaluated at param_values.

A seed derivative that reads other model symbols (d(beta_N)/d(R0_) = gamma_/N_) is only a number once the fit vector is known, so the assembly asks for the routing at the evaluated PSet. The environment is the model’s nominal parameter table overridden by the free parameters that bind a model id, then by this experiment’s condition – the same order the apply paths use (ADR-0076). A routing with no point-dependent factor returns itself, so every pre-#530 fit is untouched.

check_column_multiplicity()

Raise unless every route reads each native sensitivity column exactly once (#537).

A route’s derivative is the sum over its contributions, so a column named by two of them is added twice. route_experiment() folds same-column terms as it builds a route, which makes that impossible by construction – this re-checks it where the columns are actually summed (pybnf.gradient.assembly._raw_sensitivity_accessor(), once per experiment per evaluation), because the failure it guards against is silent: the assembled column is a clean multiple of the true derivative, so the fit walks a scaled surface and converges to a plausible wrong answer instead of erroring.

Raises PybnfError naming the free parameter, the axis and key of the doubled column, and each duplicate’s factor – the diagnostic #537 lacked when a Raia_CancerResearch2011 ic-axis column assembled at exactly 2x its central difference and did not reproduce. Never raises for a legitimate multi-path route (two condition targets seeding one species): those arrive already folded.

property is_point_dependent

Whether any chain-rule factor must be re-evaluated at the fit point (#530).

property sensitivity_ic

sensitivity_ic= for the experiment’s Simulator: every species initial-condition column any free parameter reaches, de-duplicated in declared free-parameter order.

property sensitivity_params

sensitivity_params= for the experiment’s Simulator: every parameter-axis column any free parameter reaches (a pinned = column, whose factor is a constant zero, is dropped), de-duplicated in declared free-parameter order.

class pybnf.gradient.routing.ParamRoute(free_param: str, contributions: tuple)

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

A free parameter’s derivative is the sum over its contributions – one RouteContribution per native sensitivity column it reaches. The common case is a single contribution: a parameter: free parameter bound by id (ADR-0034) reaches exactly one model column. A free parameter routed only through a condition (a per-condition estimated initial condition, ADR-0076) reaches its column through the condition target instead; and a free parameter a condition assigns to several model entities at once (a shared rate multiplier) reaches several columns, so its derivative is their sum.

.target / .key / .factor read the sole contribution of a single-column route (every bind-by-id route); reading them on a multi-column route raises – use .contributions.

classmethod single(free_param, target, key, factor)

A route with a single RouteContribution – the common bind-by-id case.

class pybnf.gradient.routing.RouteContribution(target: str, key: object, factor: float, node: tuple = None, origins: tuple = ())

One (native sensitivity column -> free parameter) term of a route.

target is PARAM (kinetic/global -> sensitivity_params), IC (species initial value -> sensitivity_ic), or NONE (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 folded into this term (#449 multiplies it into the Jacobian column); a zero factor drops it.

node carries the symbolic derivative when the factor is point-dependent – a seed whose d(entity)/d(target) reads other model symbols (#530). factor then holds its value at the routing’s build point and ExperimentRouting.at_point() refreshes it before each assembly; such a term is always requested, since a factor that merely happens to vanish at the build point must not drop the column the fit later needs.

Two contributions of one route never name the same (target, key): same-column terms are folded into one by route_experiment() (#537), so column is a route-unique identity the assembly can check against.

origins labels the chain-rule path(s) this term came from – 'bind' for the bind-by-id classification, 'ref:<target>' for a condition’s parameter-reference through that target. It carries the provenance a fold would otherwise erase: folding is what makes one-contribution-per-column structural, but it also turns two same-column terms into one term with a doubled factor, which is indistinguishable after the fact from a legitimate single term. The labels keep the two apart – one path reaching one column twice is a defect, two paths meeting on it is arithmetic – so ExperimentRouting.check_column_multiplicity() can still see through the fold. Metadata, not identity: excluded from equality/hash/repr so a routing compares by what it computes.

property column

The native sensitivity column this term reads – (target, key) (#537).

property repeated_origin

The first chain-rule path this term folded in more than once, or None.

One path reaching one native column twice is the double-count #537 is about; the fold would otherwise present it as a single term with twice the factor.

property requested

Whether this term needs its native sensitivity column computed.

class pybnf.gradient.routing.SeedTerm(target: str, key: object, node: tuple)

One (sensitivity column, d(column entity)/d(model parameter)) seeding term.

A model parameter a condition: sets may not appear in the ODE right-hand side at all: it can seed other entities’ initial values – a species initial condition (IC), or another parameter an initialAssignment derives (PARAM, beta_N = R0_*gamma_/N_). node is the symbolic d(entity)/d(parameter) tree (pybnf.gradient.derivative); it is a plain number for the common linear seed and an expression over model symbols otherwise, evaluated per fit point (#530).

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.apply_routings(model, routings)

Hand the union request over several routings to #447’s gradient path on model.

The model’s forward-sensitivity request rides the scatter and is applied at every simulate(), so it must cover every column any scored experiment reads – the union of the per-condition sensitivity_params / sensitivity_ic. (The wildtype request is not a superset once a condition routes a free parameter to a column no other experiment binds – a per-condition estimated initial condition, ADR-0076.) Capability-gated exactly as apply_routing(). Returns the applied (params, ic) lists.

pybnf.gradient.routing.classify_bound_id(name, param_ids, species_names, ic_seed_map, species_initializers=(), rhs_symbols=None, backend_ic_seeds=None)

Every sensitivity column a model id reaches: [(axis, key, node), ...], possibly empty.

The shared core of “what does this id move?”, used for both a condition target and a free parameter bound by id (ADR-0034) – they are the same question, and answering them differently is what left a bind-by-id seed with a silently-zero column (#534). An id reaches the trajectory by being a quantity the ODE reads (its own axis) and by seeding other entities’ initial values (ic_seed_map), and it may do both.

A pure initial-value seed – an id every one of whose seeds is a species IC – gets no axis of its own, provided the model confirms the ODE right-hand side never reads it (rhs_symbols, from BngsimModel.ode_rhs_symbols()). Both halves are needed. Seeding only initial conditions does not imply absence from the right-hand side: a steady-state initialAssignment seeds every species initial from the very kinetic constants that drive the rate laws, and dropping their axes deleted the entire right-hand-side half of ten derivatives in Fiedler_BMCSystBiol2016 (ADR-0097, #535). rhs_symbols is a veto on dropping, never a reason to drop.

A model that cannot say is refused (#537). ADR-0097 kept the axis there, calling the worst case one redundant sensitivity vector; that is wrong, because the axis is not zero – the backend seeds ∂x(0)/∂p into it too, so keeping it doubles the column rather than wasting a vector. With one error deleting half a derivative and the other doubling it, there is no safe default left to pick.

species_initializers is the backend-independent fallback for an id the seed map does not mention: a bare initializer (species <- p) is a unit seed, which is how classify_free_param() has always recognised an initial-condition parameter, and is what a caller that supplies no ic_seed_map at all still relies on.

Empty when the id binds nothing at all (a free sigma); the caller decides whether that is a NONE route or a refusal.

pybnf.gradient.routing.classify_condition_target(target, param_ids, species_names, ic_seed_map, rhs_symbols=None, backend_ic_seeds=None)

Classify the model entity a param-ref condition sets: return a list of (axis, key, node) seeding terms.

target is the model id a condition assignment target = free_param sets (a per-condition estimated initial condition, ADR-0076). It reaches the trajectory two ways, and a target can do both at once:

• by being a model quantity the ODE reads – a species set directly (IC on itself) or an ordinary global (PARAM on itself), factor 1; and

• by seeding other entities’ initial values (ic_seed_map) – one or more species initial conditions (I_ = I0_, S_ = N_ - I0_), and/or a parameter an initialAssignment derives from it (beta_N = R0_*gamma_/N_). Each such term carries its own d(entity)/d(target) derivative (#530).

A pure initial-value seed – a parameter that only seeds species ICs – deliberately gets no PARAM term of its own: it is absent from the ODE right-hand side, so that axis is identically zero and requesting it would only cost a sensitivity vector.

ic_seed_map maps a model parameter to its tuple of SeedTerms, or to None for a seed this routing cannot differentiate (an expression outside the arithmetic grammar, or one reaching an initial value through a rule / another assignment). That, and a target binding no sensitivity entity at all, raise GradientNotSupported – keeping a gradient/EFIM fit honest rather than emitting a silently-wrong column.

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 its parameter axis carries no separate right-hand-side path – and, because the backend seeds ∂x(0)/∂p into it, would duplicate the IC axis rather than add to it, #537). 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_values, species_initializers, condition=None, ic_seed_map=None, rhs_symbols=None, backend_ic_seeds=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_values the model’s begin parameters namespace – an {id: nominal value} mapping (any plain iterable of ids also works, but then a point-dependent seed factor has no environment to read and refuses); species_initializers the (species, initial-expr) pairs; condition the experiment’s pybnf.pset.MutationSet (None for the wildtype experiment); ic_seed_map the {model parameter -> SeedTerms} initial-value seed map (classify_condition_target()); rhs_symbols the ids the ODE right-hand side reads, or None when the model cannot say – in which case every bound id keeps its own axis (classify_bound_id()).

A parameter-reference perturbation (a per-condition estimated initial condition, ADR-0076) target = free_param composes the chain rule: the referenced free parameter reaches the trajectory through every column the condition target reaches, so it gains one RouteContribution per column – the target’s own axis, plus one per entity whose initial value the target seeds, each carrying its own d(entity)/d(target) derivative (#530). One free parameter a condition assigns to several targets at once (a shared rate multiplier) accumulates the contributions of all of them; a route’s derivative is the sum over its contributions.

A seed the arithmetic grammar cannot differentiate, a non-= parameter reference, and a target that binds no sensitivity entity all raise GradientNotSupported rather than emit a silently-wrong column.

Terms of one route that meet on the same native column are folded into a single contribution carrying their summed derivative (_fold_same_column(), #537) – the same sum the assembly would have accumulated, but leaving each column named exactly once, so a column that is named twice is unambiguously a defect and can be checked for.

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 table (id -> nominal value), (species, initial-expr) pairs, and the initial-value seed map through BngsimModel.sensitivity_entity_namespace(), plus the ids its ODE right-hand side reads through BngsimModel.ode_rhs_symbols() (the only model coupling), so the routing core stays backend-agnostic. A model exposing neither answers None to the second question, which keeps every bound id’s own axis. 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 comes from each parameter’s scale in closed form (Scale.d_inverse, priors/scale.py): 1 for LINEAR (short-circuited), ln(10)*10**u for log10, exp(u) for ln. No autodiff, so no gradient fit built on the shipped scales needs the optional pybnf[jax] extra (ADR-0087, #524); a custom scale that supplies only an inverse_jax still falls back to jax.grad (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.

The estimated-scale column generalizes past a single free parameter: noise_grad_point returns the full sum_p (dL/dp)*(dp/dtheta) vector, where each source supplies dp/dtheta (sigma_sensitivity). A single free sigma has dsigma/dname = 1 and no sim coupling, so its column is byte-identical to the historical scalar one. A composite sigma is the general chain rule dsigma/dtheta = sum_symbol (dsigma/dsymbol)*(dsymbol/dtheta): a PSet-only formula (FormulaSigma / PerMeasurementFormulaSigma, ADR-0044/0045/#505) puts dsigma/dcoeff on each coefficient’s column (the per-measurement variant binding a row token from exp_data/exp_row); a prediction-dependent sigma (sigma = sigma_abs + sigma_rel*y, a PredictionFormulaSigma, ADR-0075/0079) additionally chains dsigma/dprediction through the SAME raw_sens forward sensitivity the residual rides – so, unlike the others, it also perturbs the model-parameter columns. Every estimated scale stays on the scalar path (the retained normalizer is not a square).

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/0066). 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, which for a chain of Data-level transforms wraps the tensor accessor once per stage, each stage’s rule read in the values that stage consumed and produced, #539/ADR-0102) 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.

The two ADR-0066 primitives close the same way (#533). A floor (x' = x + rho*max(x)) is additive and separable, so ∂x'_i/∂θ = s_i + rho*s_argmax – one more term in _normalized_sensitivity. An analytic per-series scale is the one transform that is not per-point: its c* is profiled from the whole matched (sim, data) series, so the scored value is c*(θ)·s_i(θ) and its sensitivity is the product rule c*·∂s_i/∂θ + s_i·∂c*/∂θ. The series-wide ∂c*/∂θ is the closed-form derivative of the profiling condition (analytic_scale_sensitivity()) – a geometric-mean ratio for a log family, the least-squares optimum for a linear one – computed once per experiment here and shared by every point of the column. It does not vanish by the envelope theorem: the profiling criterion is family-aware but σ-unweighted, so it is not in general the objective’s own minimizer over c. Resolving which columns this experiment scales needs the experiment’s data_key, which travels as the optional 4th element of each experiments item; without it a scaled column is refused rather than silently differentiated unscaled.

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, hessian: ndarray = None)

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.

hessian: ndarray = None

The expected-Fisher / Gauss-Newton Hessian (n_param, n_param), sampling space – attached only on the EFIM trust-region path (job_type = gntr, #481/#488) by assemble_gradient_and_fisher_hessian(); None for trf / lbfgs, which never form it.

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_constraint_hessian(constraint_sets, sim_data_dict, routings, free_params)

The Gauss-Newton Hessian of the total constraint penalty w.r.t. the free parameters (layer I curvature, #481/#456), in sampling space – the constraint block the EFIM trust-region path (job_type = gntr) adds to assemble_fisher_hessian()’s data-fit Hessian, the curvature sibling of assemble_constraint_gradient().

Each constraint’s penalty is P(q(theta)) for an at-/between-time readout q, so its exact Hessian is P''(q) * outer(grad q, grad q) + P'(q) * hess q; the Gauss-Newton term drops the second-order sensitivity hess q (as the EFIM drops it for the data fit) and clamps P'' to its positive part, giving the PSD max(P''(q), 0) * outer(grad q, grad q) (penalty_curvature()). Reuses _constraint_sensitivity_accessor() for grad q exactly as the gradient does, and applies the same native -> sampling d theta/d u factor on both axes. A piecewise-linear (static hinge .con) penalty has P'' == 0, so it contributes no curvature – correct, a linear penalty has none; its pull rides the gradient. Returns zeros for an empty constraint set; raises GradientNotSupported for a penalty whose curvature this cut does not assemble (a clipped smooth penalty, an estimated constraint scale) -> the L-BFGS-B fallback.

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

Assemble the expected-Fisher / Gauss-Newton Hessian H (n_param x n_param), summed across experiments. This standalone API produces the same curvature the combined assemble_gradient_and_fisher_hessian() path feeds to the EFIM trust-region optimizer (job_type = gntr, #481/#488).

H = sum_i w_i [ kappa_i * outer(s_i, s_i)  +  sum_p I_scale_p * outer(g_i^p, g_i^p) ] where s_i = d(prediction_i)/d(theta) is the same forward sensitivity the gradient uses (through the objective’s prediction_sensitivity seam), kappa_i the per-point location Fisher (objective.location_fisher_point: (d rho/d pred)**2 for a residual-bearing Gaussian/Student-t column, the family’s location_fisher for Laplace/count), I_scale_p each estimated noise parameter’s expected Fisher, and g_i^p = d(noise_param_p)/d(theta) its scale sensitivity – the whole sum_p I_scale_p * outer(g_i^p, g_i^p) noise block returned by objective.noise_fisher_point (ADR-0080). Every rank-1 term is PSD (kappa >= 0, I_scale >= 0), so H is PSD by construction. For a bare free sigma g_i^p is the unit vector e_p (the noise parameter is model-unbound, 0 in s_i), so its block is the historical diagonal entry I_scale * outer(e_p, e_p) and the estimated-sigma Hessian is block-diagonal, exactly the Fisher predicts. For a prediction-dependent sigma (sigma = sigma_abs + sigma_rel*y, ADR-0075/0079) g_i^p also carries model-parameter columns (the scale rides the prediction), so outer(g_i^p, g_i^p) produces the genuine location↔scale coupling off the diagonal – a strict superset of the diagonal cut.

experiments is the same (sim_data, exp_data, routing[, data_key]) iterable assemble_gaussian_gradient() consumes, free_params the same ordered free- parameter list. Mirrors that assembler’s point loop exactly (same points, same raw_sens accessor, same prediction_sensitivity), so the Hessian is formed over precisely the points the gradient and objective score. The native -> sampling transform (ADR-0029) is the gradient’s d theta/d u factor applied on both axes: H <- diag(f) H diag(f). Raises GradientNotSupported for a configuration whose Fisher this cut does not assemble (a MEDIAN-centered count, a MEAN-on-log estimated scale, …), so the fit refuses the EFIM step with a pointer to the L-BFGS-B path.

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[, data_key]) items – 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. The optional 4th element is that experiment’s data_key (the suffix evaluate scores it under), needed only to resolve an analytic per-series scale (ADR-0066, #533): omitting it is byte-identical for every other fit, and refuses a scaled column rather than differentiating it unscaled. 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.

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

Assemble a GradientResult and attach its expected-Fisher Hessian in one pass.

This is the gntr objective-assembly path (#488). It produces the same residual, Jacobian, scalar gradient, and Fisher Hessian as calling assemble_gaussian_gradient() and assemble_fisher_hessian() separately, but walks each scored point only once. The shared point walk resolves the simulation row, filters missing observations, builds the raw sensitivity accessor, and calls prediction_sensitivity once before feeding both the gradient and curvature accumulators.