C — Gradient plumbing: Gaussian gradient + residual-Jacobian assembly, convention pin, native→sampling transform, FD-check

[wshlavacek]

Part of #385 (cut 1, step 3 of A→B→C). Depends on #447 and #448. Closes the common case and unblocks #386.

What

Assemble, from #447's per-Result sensitivity tensor and #448's per-experiment routing, both forms of
the gradient for the default Gaussian, LINEAR scale, fixed-σ objective, summed across experiments:

1. scalar ∂F/∂θ (for quasi-Newton / L-BFGS-B), and
2. residual vector + residual-Jacobian (for trust-region least-squares, Gradient-based local optimizers (D2D-style: TRF/Levenberg-Marquardt + L-BFGS-B) in the async Algorithm loop #386's primary path).

Convention pin

PyBNF Gaussian loss (noise/gaussian.py:52-54): data_fit = r²/(2σ²),
r = μ(pred) − additive_on.forward(obs); for default (LINEAR, MEDIAN) Gaussian, μ = pred exactly
(gaussian.py:15-17, symmetric ⇒ smooth, no trajectory centering). Therefore:

• residual ρ_i = (pred_i − obs_i)/σ_i, loss = ½‖ρ‖²
• residual-Jacobian J_ij = (1/σ_i)·(factor_j)·∂pred_i/∂θ_j (native; factor_j from B — Gradient plumbing: route free params to sensitivity_params/sensitivity_ic (bind-by-id + condition chain-rule) #448)
• ∂loss/∂θ = Jᵀρ (not 2Jᵀρ)

scipy.least_squares minimizes ½‖ρ‖² with the same ρ/J, so the residual form and the scalar
form agree by construction — the optimizer walks the surface PyBNF reports. Assembly threads through
the existing per-point objective evaluation (LikelihoodObjective.eval_point, objective.py:723) so
the gradient is the exact derivative of the loss PyBNF reports — same σ-weighting, same column selection.

Native → sampling space (once)

Apply the θ↔u Jacobian exactly once (ADR-0029): priors/scale.py already exposes a JAX-differentiable
inverse_jax per scale (scale.py:38/61/75/91), so ∂θ/∂u is one autodiff away — no hand-written
per-scale derivative, no double transform. ρ is scale-invariant; J → J·diag(∂θ/∂u).

Reaching the optimizer

The gradient must arrive through the async propose/score loop (Algorithm, algorithms/base.py:69;
make_job/got_result), not a blocking call — consumed by #386, but the assembly API here must fit
that loop (return the residual/Jacobian alongside the scored result).

Deliverables

• Objective (and Constraint stub) methods producing scalar ∂F/∂θ and residual +
  residual-Jacobian, summed across experiments, for Gaussian/fixed-σ/LINEAR.
• Native→sampling transform applied once via priors/scale.py (ADR-0029).
• FD test (acceptance gate): central differences of PyBNF's own loss(θ) vs assembled
  grad_sampling on a 2–3 parameter net model — wildtype + one conditioned experiment —
  to a tight tolerance.
• User-guide section: enabling sensitivities, the ~(1+P) cost factor, the capability gate.

Out of scope (deferred layers, each its own issue under #385)

Estimated σ (layer D); LOG10/lognormal scale objective (layer E); normalization & cumulative→incident
(layer F); asymmetric families (layer G); SBML/Antimony (layer H); constraint/comparison gradients
(layer I); pre-equilibration/steady-state (layer J); the optimizer itself (#386).

[wshlavacek]

Landed in a344665 (on main).

All deliverables met:

• Objective residual seam + Constraint stub producing scalar ∂F/∂θ and
  residual + residual-Jacobian, summed across experiments, for Gaussian/LINEAR/MEDIAN/
  fixed-σ. LikelihoodObjective.residual_point reuses the same _prediction /
  _noise_values seams as eval_point, so ½‖ρ‖² == evaluate exactly; the base
  objective and Constraint.penalty_gradient (layer I) stub-raise the new
  GradientNotSupported. assemble_gaussian_gradient returns
  {residual ρ, jacobian J, gradient = Jᵀρ} (the pin: ∂loss/∂θ = Jᵀρ, not 2Jᵀρ).
• Capability gate: Gaussian × LINEAR × MEDIAN × fixed-σ, no measurement layer /
  per-measurement / cumulative — every other config raises a pointed
  GradientNotSupported naming its deferred layer (D–J).
• Native→sampling transform applied once (ADR-0029): J[:,j] *= dθ/du, ρ invariant;
  LINEAR short-circuits (no jax), log scales autodiff inverse_jax via a lazily-imported
  pybnf[jax].
• FD acceptance gate: tests/test_gradient_assembly.py (15 tests) — central
  differences of PyBNF's own loss vs the assembled gradient on the decay net (k
  parameter-axis + S0 initial-condition-axis, wildtype + a k*4 condition, linear &
  log10) matching to ~1e-7. Added to the pre-push bngsim gate.
• User-guide section: docs/gradient_fitting.rst (enabling sensitivities, the
  ~(1+P) cost factor, the capability gate).

Out-of-scope coverage layers D–J remain open under #385. The FD test surfaced a
pre-existing bngsim-backend bug — a free parameter that only seeds a species initial
condition is a silent no-op via set_param — spun off as #450 (it blocks #386's IC fits).