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.

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 the target’s own sensitivity column (classify_condition_target()): IC when target seeds a species initial value (chain-rule factor d(IC)/d(target) = 1 for a bare initialAssignment), PARAM for an ordinary global (e.g. a shared rate multiplier). Because one free parameter may be assigned to several targets in one condition, a route is a sum over its contributions (ParamRoute).

Cut-1 scope

IC routing matches a bare initializer (species <- p directly, or target = p where target bares a species IC). A parameter that both appears in the ODE RHS and seeds an initial value, or a species seeded by a non-trivial expression (2*p), is a later layer – classify_condition_target() raises GradientNotSupported for the non-bare seed rather than emitting a parameter-dependent factor, keeping the clean param/ic split honest.

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: every species initial-condition column any free parameter reaches with a non-zero factor, 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 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, 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)

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.

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_condition_target(target, param_ids, species_names, ic_seed_map)

Classify the model entity a param-ref condition sets: return (axis, key, factor).

target is the model id a condition assignment target = free_param sets (a per-condition estimated initial condition, ADR-0076). A model parameter that seeds a species’ initial value routes to the IC axis keyed by that species (checked first – an IC seed is also a begin parameters id, but only its IC axis is non-zero); a species set directly routes to IC on itself; any other model parameter routes to PARAM. factor is d(model entity)/d(model parameter) – 1 for a bare seed or a direct set.

ic_seed_map maps a model parameter to the species whose initial value it bares (a bare initialAssignment species = <param>); the value None marks a parameter that seeds a species IC through a non-bare expression (a parameter-dependent d(IC)/d(param) this bind-by-id routing does not model). Raises GradientNotSupported for that non-bare seed and for a target that binds no sensitivity entity at all – keeping a gradient/EFIM fit honest rather than emitting a silently-zero 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 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, ic_seed_map=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); ic_seed_map the {model parameter -> species} bare-initialAssignment map (classify_condition_target()).

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 the condition target’s own sensitivity column, so it gains a RouteContribution on that column (IC when the target seeds a species IC, PARAM for an ordinary global). One free parameter a condition assigns to several targets at once (a shared rate multiplier) accumulates one contribution per target – its derivative is their sum. A target whose d(IC)/d(param) is not a bare 1 (a non-bare initialAssignment), or that binds no sensitivity entity, raises GradientNotSupported rather than emitting a silently-zero column.

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, (species, initial-expr) pairs, and the bare-initialAssignment seed map 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.

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). 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, 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 (fit_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 (fit_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 (fit_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) 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) 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.

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.