Pluggable optimizer protocol with scipy wrappers for gradient-based methods

[wshlavacek]

Summary

Add a pluggable optimizer protocol to PyBioNetFit, with adapter wrappers
for scipy.optimize.minimize and scipy.optimize.least_squares. The
current optimizer suite is metaheuristic only (genetic algorithm,
particle swarm, simulated annealing, scatter search, etc.); gradient-based
methods are not currently available.

Depends on #385.

Current state

• pybnf/algorithms.py contains hand-written implementations of the
  metaheuristic optimizers. There is no shared optimizer protocol or
  registry; each algorithm is its own class with its own conventions.
• All current algorithms are derivative-free; none consume the gradient
  surfaced by Gradient plumbing: surface objective gradients via the BNGsim backend #385.
• The PyBioNetFit configuration schema selects an algorithm by name via
  the fit_type keyword; supported values are the existing metaheuristic
  names.

Scope of this issue

Add three pieces to PyBioNetFit:

1. Optimizer protocol and registry under pybnf/optimizers/. Define
   a common Optimizer base class and an OptimizationResult return
   type; provide a name → class registry so new optimizers can be added
   without modifying the dispatch code.
2. scipy.optimize.minimize wrapper supporting the scalar-objective
   methods: L-BFGS-B, BFGS, trust-constr, Newton-CG, trust-ncg,
   SLSQP, COBYLA, Nelder-Mead, Powell, CG. A single wrapper
   class parameterized by the method string.
3. scipy.optimize.least_squares wrapper supporting the residual-
   shape methods: trf, dogbox, lm. A separate wrapper class
   because the objective interface differs (residual vector rather than
   scalar).

Out of scope

• Third-party optimizer libraries (fides, cyipopt for IPOPT,
  nlopt, etc.). Each is a separate follow-up.
• Algorithms with hand-rolled implementations specific to PyBioNetFit;
  this issue covers wrappers around external libraries only.
• Hyperparameter-tuning infrastructure for optimizer options.

Proposed interfaces

Core protocol

# pybnf/optimizers/base.py

@dataclass
class OptimizationResult:
    x: np.ndarray          # best parameter vector
    fun: float             # best objective value
    success: bool
    n_fev: int             # objective evaluations
    n_jev: int = 0         # gradient evaluations (0 for derivative-free)
    message: str = ""
    history: list | None = None   # optional per-iteration log

class Optimizer(ABC):
    """Common interface for local optimizers wrapped from external libraries."""

    @abstractmethod
    def fit(self,
            objective: Callable[[np.ndarray], float],
            x0: np.ndarray,
            bounds: Sequence[tuple[float, float]],
            *,
            gradient: Callable[[np.ndarray], np.ndarray] | None = None,
            max_iter: int = 1000,
            tol: float = 1e-6,
            callback: Callable[[np.ndarray], None] | None = None,
            ) -> OptimizationResult:
        ...

For residual-shape optimizers (least_squares family), an alternative
signature is exposed as fit_residuals taking a residual callable
returning a vector, with an optional Jacobian callable returning a
matrix.

Registry

# pybnf/optimizers/registry.py

_REGISTRY: dict[str, type[Optimizer]] = {}

def register(name: str):
    def deco(cls): _REGISTRY[name] = cls; return cls
    return deco

def get_optimizer(name: str, **kwargs) -> Optimizer:
    if name not in _REGISTRY:
        raise ValueError(f"Unknown optimizer: {name!r}; "
                         f"available: {sorted(_REGISTRY)}")
    return _REGISTRY[name](**kwargs)

scipy.optimize.minimize wrapper

# pybnf/optimizers/scipy_min.py

@register("scipy_minimize")
class ScipyMinimizeOptimizer(Optimizer):
    SUPPORTED = frozenset({
        "L-BFGS-B", "BFGS", "trust-constr", "Newton-CG", "trust-ncg",
        "SLSQP", "COBYLA", "Nelder-Mead", "Powell", "CG",
    })

    def __init__(self, method: str = "L-BFGS-B", **options):
        if method not in self.SUPPORTED:
            raise ValueError(f"method {method!r} not supported")
        self.method = method
        self.options = options

    def fit(self, objective, x0, bounds, *,
            gradient=None, max_iter=1000, tol=1e-6, callback=None):
        from scipy.optimize import minimize
        res = minimize(
            fun=objective, x0=np.asarray(x0), jac=gradient,
            method=self.method, bounds=bounds,
            options={"maxiter": max_iter, "gtol": tol, **self.options},
            callback=callback,
        )
        return OptimizationResult(
            x=res.x, fun=float(res.fun), success=bool(res.success),
            n_fev=int(res.nfev),
            n_jev=int(getattr(res, "njev", 0)),
            message=str(res.message),
        )

scipy.optimize.least_squares wrapper

# pybnf/optimizers/scipy_lsq.py

@register("scipy_least_squares")
class ScipyLeastSquaresOptimizer(Optimizer):
    SUPPORTED = frozenset({"trf", "dogbox", "lm"})

    def __init__(self, method: str = "trf", **options):
        if method not in self.SUPPORTED:
            raise ValueError(f"method {method!r} not supported")
        self.method = method
        self.options = options

    def fit_residuals(self, residuals, x0, bounds, *,
                      jacobian=None, max_iter=1000, tol=1e-6):
        from scipy.optimize import least_squares
        res = least_squares(
            fun=residuals, x0=np.asarray(x0), jac=jacobian or "2-point",
            method=self.method, bounds=bounds,
            max_nfev=max_iter, gtol=tol, **self.options,
        )
        return OptimizationResult(
            x=res.x, fun=float(res.cost), success=bool(res.success),
            n_fev=int(res.nfev),
            n_jev=int(getattr(res, "njev", 0)),
            message=str(res.message),
        )

Configuration integration

In a .conf file:

fit_type = scipy_minimize
optimizer_method = L-BFGS-B
optimizer_options = maxiter=5000, gtol=1e-8

The existing fit_type keyword is extended to accept any name in the
optimizer registry. When fit_type resolves to a registered optimizer
class, the algorithm dispatch in algorithms.py constructs the
optimizer via get_optimizer(fit_type, method=optimizer_method, **optimizer_options) and drives it with the objective callable and the
gradient callable (when available from #385).

Edge cases and design considerations

• No gradient available. If a scalar-objective method that benefits
  from a gradient (L-BFGS-B, BFGS, trust-constr, Newton-CG,
  etc.) is selected and no gradient is supplied, scipy falls back to
  finite differences. This works but is expensive. Document the
  behavior; do not silently ignore the cost.
• Method does not support bounds. Some scipy methods (Newton-CG,
  CG, BFGS) do not accept bounds. Raise a clear error rather than
  silently dropping bounds.
• Method requires a Hessian. trust-ncg, trust-exact,
  Newton-CG accept a Hessian callable. The current scope does not
  surface Hessians (see Gradient plumbing: surface objective gradients via the BNGsim backend #385's out-of-scope statement on second-order
  sensitivities); these methods are still usable via scipy's
  Hessian-by-finite-differences fallback. Document the cost.
• scipy.optimize.least_squares Jacobian. The residual-shape
  wrapper accepts a Jacobian callable returning the full Jacobian
  matrix (n_residuals × n_params). The residual decomposition for the
  combined SSR + qualitative-loss objective is a follow-up; the wrapper
  itself just adapts the interface.
• Multistart orchestration. Out of scope here. The existing
  multistart logic in algorithms.py continues to drive starting
  points; each restart calls the wrapped optimizer once.
• Optimizer-level callbacks. The standardized callback(x)
  interface is the lowest-common-denominator across scipy methods. Some
  methods provide richer callback interfaces (e.g., trust-constr's
  OptimizeResult-style callback); these can be surfaced in later
  refinements if needed.

Deliverables

• pybnf/optimizers/base.py: Optimizer ABC + OptimizationResult
  dataclass.
• pybnf/optimizers/registry.py: name → class registry with
  register decorator and get_optimizer lookup.
• pybnf/optimizers/scipy_min.py: ScipyMinimizeOptimizer wrapping
  scipy.optimize.minimize with the 10 supported methods.
• pybnf/optimizers/scipy_lsq.py: ScipyLeastSquaresOptimizer
  wrapping scipy.optimize.least_squares with trf, dogbox, lm.
• Integration in algorithms.py: fit_type dispatch consults the
  optimizer registry when a registered name is supplied.
• Configuration-schema support for optimizer_method and
  optimizer_options keywords.
• Unit tests: each supported scipy method on a small smooth
  objective with known minimum; verify convergence and that
  OptimizationResult fields are populated.
• Documentation: a "Gradient-based optimizers" section in the user
  guide listing the supported methods and configuration syntax.
• No behavioral change for existing metaheuristic algorithms.

Backward compatibility

Purely additive. The existing fit_type values for hand-rolled
algorithms continue to resolve to the existing algorithm classes. New
fit_type values for registered optimizers route through the new
adapter layer. No existing API method or configuration keyword changes
meaning.