aquakin.calibrate#

aquakin.calibrate(reactor, C0, observations, t_obs, free_params, *, transforms=None, initial_params=None, observed_species=None, time_unit=None, loss='mse', sigma=None, priors=None, use_priors=True, diff=DifferentiationConfig(mode='reverse', method='through_solve', check_finite=True, adjoint_max_steps=100000, adjoint_low_memory=False), optimizer=OptimizerConfig(method='lbfgsb', n_starts=1, jitter=0.5, jitter_schedule=None, seed=0, max_iter=500, tol=1e-06, param_halfwidth=None), laplace=True, free_ic=None, _compiled_cache=None)[source]#

MAP fit with optional Laplace posterior approximation.

Parameters:

  • reactor (Reactor) – Batch reactor (or anything with a compatible solve). Same usage contract as fit().

  • C0 (jnp.ndarray or list of jnp.ndarray) – Initial concentration vector. Pass a list of vectors for a joint multi-batch fit: each entry is one dataset’s initial state, the batches share the parameter vector and prior, and their data terms are summed. observations and t_obs must then be matching lists.

  • observations (jnp.ndarray or list of jnp.ndarray) – Observed values, shape (n_t,) for a single species or (n_t, n_observed). In multi-batch mode, a list of such arrays (one per dataset; the datasets may have different n_t).

  • t_obs (jnp.ndarray or list of jnp.ndarray) – Observation times, shape (n_t,). C0 is taken at t=0; the solver integrates from 0 to t_obs[-1]. In multi-batch mode, a list of time grids, one per dataset. In the model’s native time unit unless time_unit is given.

  • free_params (list[str]) – Namespaced parameter names to calibrate. Others held fixed.

  • transforms (dict[str, str], optional) – Per-parameter transform. Keys may be any subset of free_params; unspecified entries fall back to the parameter’s declared transform on the model (default "none").

  • initial_params (jnp.ndarray, optional) – Starting parameter vector. Defaults to reactor.model.default_parameters().

  • observed_species (list[str], optional) – Species names corresponding to columns of observations. If None, every model species is taken to be observed.

  • time_unit (str, optional) – The time unit t_obs is expressed in ("s", "min", "h", "d"), matching BatchReactor.solve(). Every dataset’s t_obs is converted into the model’s native (rate-constant) time unit before the solve, so an hour-valued t_obs carried over from a solve(time_unit="h") run is interpreted correctly rather than as native-unit days (the silent 24x time-axis compression this guards against). Default None interprets t_obs in the native unit. The fitted rate constants are always in native units.

  • loss ({"mse", "wmse", "nll"}, optional) – Loss function.

  • sigma (jnp.ndarray, optional) – Per-observation standard deviation for "wmse" / "nll". Scalar, (n_observed,), or (n_t, n_observed). In multi-batch mode, either a single value/array shared across datasets or a list with one entry per dataset.

  • priors (dict[str, tuple[float, float]], optional) – Gaussian priors as name -> (mean, std) in physical space, added to the objective as 0.5 * sum(((p - mean) / std) ** 2). Overrides any prior declared on the model for the same parameter. Only entries whose name is in free_params are used.

  • use_priors (bool, optional) – If True (default), parameters whose model declaration carries a prior: block contribute their Gaussian prior to the objective (for the free parameters), in addition to any passed via priors. Set False to ignore the model-declared priors. Priors regularise otherwise non-identifiable parameter combinations toward literature values; for a proper Bayesian MAP / posterior, combine them with loss="nll" and a measurement sigma so the data term is a true negative log-likelihood (the prior curvature then enters the Laplace covariance automatically).

  • optimizer (OptimizerConfig, optional) – Optimiser backend and multistart settings; see OptimizerConfig for the fields (method {“lbfgsb”, “gauss_newton”}, n_starts, jitter / jitter_schedule, seed, max_iter, tol, param_halfwidth). "gauss_newton" minimises the residual vector with SciPy least_squares (its Jacobian by AD, direction from diff.mode) and is markedly more robust on the multimodal landscapes of stiff reaction-model fits; n_starts > 1 is a deterministic, seed-reproducible multistart that keeps the lowest-loss optimum. Default OptimizerConfig() (a single L-BFGS-B start).

  • laplace (bool or LaplaceConfig, optional) – Laplace covariance approximation at the MAP. True (default) computes it with default settings, False disables it, and a LaplaceConfig computes it with tuning (method {“fd”, “gauss_newton”}, ridge, eig_keep, fd_step, dtmax). It is interpretable as a Bayesian posterior only when loss="nll" with a calibrated sigma (the loss then IS a proper Gaussian negative log-likelihood); for "mse" / "wmse" it is the inverse loss curvature – the right shape but not the right absolute scale. See fit() if you only need point estimates.

  • free_ic (FreeICConfig, optional) – If given, also fit the initial concentration of the named species (per dataset, in log space, box-bounded by FreeICConfig.bounds, with an optional log-space prior via prior_log_std) while the rate parameters stay shared – useful when an unmeasured initial state (e.g. a biofilm reservoir) is not known. The fitted pools are returned on CalibrationResult.C0_fitted / ic_named; the Laplace posterior is taken over the rate parameters with the pools held at their optimum. Default None (no free ICs).

  • diff (DifferentiationConfig, optional) – How the data-term gradient / residual Jacobian is formed. mode ({“reverse”, “forward”}) is the autodiff direction; method ({“stable”, “through_solve”}) is how the reverse adjoint is formed. mode="reverse", method="through_solve" (the calibrate default) differentiates through the diffrax solve (RecursiveCheckpointAdjoint); for a stiff model this reverse pass goes non-finite above a step-size threshold, so the reactor must carry a dtmax cap. mode="reverse", method="stable" replaces the integrator’s adjoint with an explicit per-step transposed-stage solve (esdirk_adjoint_solve()) – a robust adaptive ESDIRK forward whose backward is finite with no cap (batch reactors only; its gradient agrees with the capped through_solve path to the optimiser tolerance). mode="forward" uses jacfwd and builds a forward-capable reactor adjoint internally (forward-mode AD stays finite at any step – the fix for a stiff model whose reverse adjoint overflows); pair it with optimizer="gauss_newton" and method="through_solve" (forward + method="stable" is rejected, the stable adjoint being reverse-only). check_finite (default True) raises a friendly error if the start-point gradient is non-finite. adjoint_max_steps is the (allocated) forward step count for the method="stable" solve – the backward scan walks this whole saved-trajectory buffer, so set it to a tight upper bound on the step count (ignored for through_solve). adjoint_low_memory recomputes each step’s stages in the backward pass instead of saving the ~n_stages × dense-stage buffer – memory for compute, with the gradient unchanged.

  • _compiled_cache (dict | None)

Return type:

CalibrationResult