steady_state() never receives the compiled RHS, and finite-differences all of dY_ss/dp — including the Jacobian the model already has analytically

[wshlavacek]

Simulator.steady_state() and steady_state_batch() run the interpreted RHS
unconditionally, and their sensitivity matrix is built entirely from finite
differences. Both are invisible from Python: nothing warns, and the result comes
back looking like every other sensitivity result.

Evidence

1. The codegen field exists, is exposed, and is read by nobody.

SteadyStateOptions carries std::string codegen_so_path (include/bngsim/types.hpp:782),
and pybind exposes it. But:

• Simulator.steady_state() builds a SteadyStateOptions and sets tol,
  max_time, method, rtol, atol, max_steps, jacobian and
  sensitivity_params — never codegen_so_path, even when the Simulator holds a
  .so in self._codegen_so_path. steady_state_batch()'s _run_one does the
  same.
• src/steady_state.cpp contains exactly one occurrence of the string "codegen",
  a comment at line 947. The option is never read on the C++ side either.

So a Simulator built with codegen=True solves for steady state on the same
interpreted path as one built with codegen=False. There is also no
codegen_c_source field at all, so the MIR JIT backend is not even representable
here.

2. dY_ss/dp is finite-differenced twice over.

compute_ss_sensitivity (src/steady_state.cpp:703) forms dY_ss/dp = -J⁻¹·(∂f/∂p) as:

• J by finite differences — n_species + 1 interpreted RHS evaluations, one per
  column, with the comment "We use finite differences (works for all rate law
  types)" (lines 731–753). The model's complete analytical Jacobian is not
  consulted, even though the newton path's KINSOL polish in this same file uses
  it, and jacobian="auto" selects it everywhere else.
• ∂f/∂p by finite differences — one more RHS evaluation per parameter
  (lines 755–781), by mutating the parameter in place.

Both use a fixed sqrt(eps) step. On a 1281-species model that is ~1300
interpreted RHS evaluations to assemble one Jacobian, and the codegen analytical
∂f/∂p — the same derivative the CVODES sensitivity RHS emits — is not used.

Why it matters

GH #214 retired the interpreted finite-difference sensitivity path for the
time-course solve on the grounds that its ~sqrt(eps) noise cannot support tight
tolerances, and made Simulator(..., sensitivity_params=...) refuse rather than
degrade. Issue #62 is the same defect surviving in compute_all_sensitivities.
The steady-state path is the third instance, and the only one where the
finite-difference derivative is the design, not a fallback: it will not be
fixed by the #62 change.

The docs are half-right about this. docs/user-guide/steady-state.md:142 says
∂f/∂p is computed by finite differences, but describes J as "the Jacobian at
steady state", which reads as the analytical one the rest of the library uses.

Suggested fix

Three separable pieces, in increasing order of work:

1. Pass self._codegen_so_path through in steady_state() and
   steady_state_batch(), and honor it in steady_state.cpp (and add a
   codegen_c_source field for the JIT backend). This alone speeds the state
   solve up and gives the sensitivity code a compiled RHS to difference.
2. Use the analytical Jacobian for J when analytical_jacobian_complete is set,
   falling back to finite differences otherwise — the same rule jacobian="auto"
   already applies, and the same matrix the newton polish already uses.
3. Use the codegen analytical ∂f/∂p when available, which is what makes the
   result consistent with the GH #214 policy for the time-course path.

Until (2) and (3) land, the docs should say plainly that both factors are finite
differences.