Gradient plumbing: surface objective gradients via the BNGsim backend

[wshlavacek]

Summary

Add the infrastructure required to support gradient-based optimization in
PyBioNetFit. The current optimizer suite is metaheuristic only (genetic
algorithm, particle swarm, simulated annealing, scatter search, etc.);
gradient-based methods are not currently available.

This issue covers the gradient-plumbing prerequisite: surfacing the
gradient of the fitting objective with respect to free parameters from
the simulation backend through to the optimizer layer. A follow-up issue
will cover optimizer wrappers themselves.

Current state

• pybnf/bngsim_model.py calls bngsim.Simulator(...) for simulation
  and returns simulation trajectories. It does not request
  sensitivity_param=[...] or sensitivity_ic=[...] from BNGsim, and
  does not expose any per-condition Jacobian information to its caller.
• pybnf/objective.py computes scalar objective values from simulation
  trajectories. It has no gradient method.
• pybnf/constraint.py computes per-constraint scalar penalties via
  get_static_penalty and per-constraint log-likelihoods via
  get_log_likelihood. Neither has a companion gradient method.
• pybnf/algorithms.py has no pathway for passing gradient information
  to an optimizer.

BNGsim itself supplies forward sensitivities (∂g/∂θ_p along the
trajectory of each output variable g) via the sensitivity_param=[...]
and sensitivity_ic=[...] kwargs on its Simulator constructor, with
results available through the Result.sensitivities and
Result.xr.sensitivities accessors. The plumbing work is on the
PyBioNetFit side.

Scope of this issue

Add three pieces to PyBioNetFit:

1. Model-layer Jacobian surfacing. Extend BngsimModel (and the
   SBML/Antimony variants in bngsim_sbml_model.py,
   bngsim_antimony_model.py) so that when a list of fitting parameter
   names is supplied at simulation time, the returned simulation-data
   object includes per-time-point Jacobian entries ∂g/∂θ_p for each
   requested output variable.
2. Objective-layer gradient assembly. Add gradient methods to
   objective.py and constraint.py that consume the model-layer
   Jacobians and assemble ∂F/∂θ for the combined objective
   F = F_quant + F_qual.
3. Algorithm-layer dispatch. Provide a uniform way for an algorithm
   to request the gradient at a parameter point and receive it. Existing
   metaheuristic algorithms ignore this; the gradient pathway is
   exclusively for gradient-based algorithms that will be added in
   follow-up work.

Out of scope

• Specific optimizer implementations or wrappers.
• Second-order sensitivities (Hessians). The forward-sensitivity
  machinery in BNGsim is first-order; this issue does not introduce
  second-order plumbing.
• Adjoint-mode sensitivities. Forward-mode only.

Proposed interfaces

Model layer

class BngsimModel(Model):
    def execute(self, params: dict, *,
                sensitivity_params: list[str] | None = None,
                sensitivity_ic_params: list[str] | None = None) -> SimulationData:
        ...

When sensitivity_params or sensitivity_ic_params is supplied, the
returned SimulationData exposes a sensitivities accessor:

data.sensitivities[param_name]  # indexed by output variable and timepoint

If neither is supplied, behavior is unchanged from current.

Objective layer

class Objective:
    def evaluate_gradient(self, sim_data, params: dict) -> dict[str, float]:
        """Return ∂F/∂θ_p for each parameter in params."""
        ...

The implementation composes per-data-point sensitivities into the
gradient of the SSR + qualitative-loss combination, summing
contributions across all data points and constraints.

Constraint layer

class Constraint:
    def get_static_penalty_gradient(self, sim_data_dict, ...) -> dict[str, float]: ...
    def get_log_likelihood_gradient(self, sim_data_dict, ...) -> dict[str, float]: ...

Each returns the gradient of the corresponding scalar method with
respect to the requested fitting parameters, via the chain rule from the
model-layer Jacobians.

Algorithm layer

Algorithms that want gradients call:

sim_data = model.execute(params, sensitivity_params=self.gradient_params)
loss = objective.evaluate(sim_data, params)
grad = objective.evaluate_gradient(sim_data, params)

Algorithms that do not (existing metaheuristics) call execute()
without sensitivity kwargs and never call evaluate_gradient(). No
change to existing behavior.

Edge cases and design considerations

• Parameters affecting only initial conditions are handled via
  BNGsim's sensitivity_ic=[...] rather than sensitivity_param=[...].
  The plumbing should accept both kinds and route each parameter name to
  the appropriate BNGsim kwarg. A user-facing way to declare which
  parameters are initial-condition parameters may be needed; existing
  model definitions may already distinguish these.
• Log-transformed parameters. PyBioNetFit fits in log space for many
  parameter types. BNGsim returns Jacobians in the model's native
  parameter space; for log-space optimization, multiply by ∂θ/∂(ln θ) = θ
  per the chain rule.
• Multiple-condition fits. Each experimental condition becomes a
  separate simulation; per-condition Jacobians compose into the total
  gradient additively.
• Comparison-difference constraints with δ_i = g(a_i, θ) − g(b_i, θ)
  require taking the difference of two per-condition Jacobians:
  ∇_θ δ_i = ∇_θ g(a_i, θ) − ∇_θ g(b_i, θ).
• Cost of sensitivities. Forward sensitivities increase simulation
  cost by approximately a factor of (1 + P), where P is the number of
  fitting parameters. This is the expected cost for users opting into
  gradient-based optimization.

Deliverables

• BngsimModel.execute() accepts sensitivity_params and
  sensitivity_ic_params kwargs and returns sensitivities when
  requested.
• Equivalent updates to BngsimSbmlModel.execute() and
  BngsimAntimonyModel.execute().
• Objective.evaluate_gradient() method.
• Constraint.get_static_penalty_gradient() and
  Constraint.get_log_likelihood_gradient() methods.
• Unit tests verifying analytical gradients against a
  finite-difference reference on a small example.
• Documentation: a section in the user guide explaining how to
  enable sensitivities and what cost to expect.
• No behavioral change for existing metaheuristic-based fits.

Backward compatibility

This issue is purely additive. No existing API method changes signature;
the new sensitivity machinery is opt-in via new keyword arguments.
Existing fits using metaheuristic algorithms continue to work unchanged.