Backend Design

This page is the contributor map for ScientificFitting's numerical core. User-facing signatures and defaults are in the API Reference; the statistical derivations are in Mathematics and Statistics.

The dependency direction is one-way:

\[\begin{aligned} \text{public constructor} &\longrightarrow \text{validated problem} \\ &\longrightarrow \text{objective and compatible solver} \\ &\longrightarrow \text{result} \\ &\longrightarrow \text{reports, diagnostics, profiles, and plots}. \end{aligned}\]

Later stages may inspect an earlier result, never reconstruct or mutate the statistical problem behind it.

Two Problem Families

A residual vector and a general likelihood do not expose the same information, so there are two normalized problem types.

Both families share parameter bounds, fixed parameters, Gaussian parameter terms, nonlinear constraints, multistart selection, local covariance, diagnostics, and profile refits. They stay separate where their data contracts differ: a generic likelihood need not have x-y residuals, predictions, or a fit curve.

One Fit, Step By Step

1. Normalize and validate the scientific input

Problem construction and the public fit entry points copy numeric inputs into stable storage and reject mismatched dimensions, non-finite observations or starting values, non-positive standard deviations, invalid covariance matrices, inconsistent parameter indices, and fixed values outside declared bounds before solver dispatch.

2. Map full parameters to optimizer coordinates

The scientific model always sees the complete parameter vector; fixed parameters are removed only from the optimizer-visible vector:

\[q_{\mathrm{free}} \xrightarrow{\text{expand}} p_{\mathrm{full}}.\]

Bounds are reduced to the same free coordinates, and nonlinear constraint callbacks still receive p_full. After fitting, free covariance and Jacobian blocks are embedded back into full parameter order. Ordinary fits, multistart candidates, profiles, and contours reuse this one mapping.

3. Prepare reusable evaluation state

Static Gaussian uncertainty is prepared outside repeated scalar objective evaluations:

• diagonal errors become inverse standard deviations and a log determinant,
• dense covariance becomes a Cholesky factor and a log determinant,
• sparse covariance keeps its sparse Cholesky factor,
• WhiteningOperator keeps the supplied matrix-free operation and determinant,
• correlated parameter constraints are factorized once.

This state lives in FitEvaluationCache (Gaussian) and LikelihoodEvaluationCache (likelihood); the LsqFit path prepares equivalent static weights for its native residual interface. Parameter-dependent covariance — effective x errors and model-relative components — is recomputed at each parameter point, not cached.

4. Construct exactly one objective

For a Gaussian problem, cost=:auto selects static chi-square when the covariance is parameter independent:

\[\chi^2(p)=r(p)^\mathsf{T}V^{-1}r(p).\]

Here $r(p) = y - f(x, p)$ is the residual vector of the model $f$ at parameters $p$, $V$ is the observation covariance matrix, and $n$ is the number of observations; $V(p)$ below denotes the same covariance when it depends on the fitted parameters.

When the effective covariance depends on the fitted parameters, it selects the normalized Gaussian objective:

\[-2\log L(p) =n\log(2\pi)+\log\det V(p)+r(p)^\mathsf{T}V(p)^{-1}r(p).\]

The determinant term then changes the optimum and cannot be dropped; see Full Gaussian Likelihood. Likelihood problems provide their data objective directly on the $-2\log L$ scale. The shared parameter layer adds Gaussian priors and correlated parameter constraints; bounds and fixed parameters restrict the parameter space, never as hidden penalty terms.

5. Dispatch only to a compatible solver

Derivative selection is stored in the problem. derivatives=:auto keeps forward-mode automatic differentiation (ForwardDiff); derivatives=:finite uses central finite differences for foreign callbacks that accept ordinary floats but not ForwardDiff duals, and applies to optimization, constraint derivatives, post-fit covariance, predictions, and profile/contour refits. Explicit model Jacobians and x-derivatives take precedence. Tolerance and budget defaults are tabulated under Solver Control: the finite-mode 1e-6 versus 1e-10 tol split and the per-solver meaning of maxiters. NativeMinuit's EDM default of 0.1 is documented under Solver Adapters. On the LsqFit path, maxiters maps to maxIter and the single tol to both x_tol and g_tol. Explicit tol values are preserved, and the result keeps the solver's actual convergence flag, so an iteration limit never implies success.

Finite mode differentiates the whole objective, including any parameter-dependent covariance and its log determinant; differentiated models must be evaluable in a neighborhood of each evaluation point, including near declared bounds. Pointwise x-error propagation uses four vectorized model calls for a central fourth-order estimate of every df/dx with step $h_i = \epsilon^{1/5}\max(|x_i|,1)$, with machine epsilon $\epsilon = \mathrm{eps}(\mathrm{Float64}) \approx 2.2\times10^{-16}$; analytic x_derivative callbacks bypass it.

The Python bridge converts inputs once and enters the fit through invokelatest, so Julia gives up compile-time specialization exactly once per fit rather than inside the numerical loop. Foreign callbacks are wrapped in a private _TypedCallback{R} whose declared result type R lets the solver code stay compiled, at the price of one dynamic call per model evaluation; different Python models therefore share Julia solver specializations, while native Julia models bypass the wrapper. Models and x-derivatives return vectors, allocating parameter Jacobians return matrices. A PrecompileTools workload caches Gaussian/Poisson kernels, prediction, and reporting without starting Python or loading plotting.

Density helpers accept vectorized=true without changing their scalar default: one batched density call over all observations per objective evaluation instead of one scalar call per point, QuadGK's BatchIntegrand with per-bin adaptive error control, dual-number-preserving batch buffers reused across bins, and the same density closure in profile refits. An analytic bin integral supplied through fit_histogram_model avoids quadrature entirely.

Solver selection follows the represented problem rather than a speed preference:

There is no backend forcing: features outside static least squares route to the scalar solver automatically. Backend selection may change how the objective is minimized, never which objective is minimized.

The Solver Extension Boundary

Both problem families enter src/solvers.jl for scalar minimization; the specialized LsqFit residual path remains in src/fit.jl. An explicit solver::AbstractFitSolver overrides automatic numerical selection, not the statistical model. An adapter implements:

• solver_capabilities: support for bounds/constraints and required derivatives.
• solve_fit: minimize the prepared OptimizationProblem and return FitSolverResult in its free coordinates, with native result and actual status.
• Optionally default_fit_tolerance: supply the solver's stopping tolerance when the user omits tol; the core stores it before multistart or profiling.

The core owns objective construction, parameter mapping, capability checks, derivative policy, multistart ranking, statistical summaries, and profile refits. Adapters do not add priors, reinterpret error scales, or silently drop constraints; unknown iteration counts stay missing. Solver settings persist when a profile changes the number of free parameters, so NLopt adapters take algorithm enums rather than dimension-bound native Opt instances.

OptimizationSolver uses SciML's algorithm traits and termination codes. ScientificFittingNativeMinuitExt loads only with NativeMinuit and supplies the same cost and derivative policy to MIGRAD with errordef=1; its native result keeps covariance and failure evidence behind a non-specializing result field, so one result type serves all models and solvers while the numerical kernels still specialize. Conventions: Solver Adapters.

The likelihood solver and parameter_covariance keywords are independent options stored in FitOptions and retained by profile refits; their contracts, including Nelder-Mead's evaluation budget and Hessian-free default, are in Minimization And Local Errors.

One multistart ranking rule serves both families: the candidate with the lowest finite cost wins, and convergence status only breaks exact cost ties. The cost is the fitted quantity, so a formally converged stop at a worse minimum never displaces a better minimum whose solver stopped without a convergence code; a non-converged winner is returned with converged=false instead of being silently replaced. If no candidate produces a finite result, the last error is raised.

ForwardDiff tags must not be shared between precompiled fits and unseen models with nested x derivatives (ForwardDiff #714); the startup gate checks exactly this in a fresh process. CHOLMOD's sparse solves reject ForwardDiff duals, so static sparse covariance requires derivatives=:finite with the general optimizer or an AD-compatible WhiteningOperator; the least-squares path needs no override, Python uses finite derivatives consistently (including sparse fits with bounds and profile refits), and sparse covariance components are validated through their stored entries, without a dense copy.

6. Build the result once

FitResult and LikelihoodFitResult are the numerical source of truth: selected minimum, solver status, parameter estimates, local covariance, correlations, statistics, and diagnostics. Gaussian x-y results additionally retain predictions, raw and weighted residuals, and the weighted Jacobian.

For static least squares, local covariance comes from the weighted Jacobian; otherwise from the objective Hessian on the $-2\log L$ scale, evaluated once and reused for diagnostics. With parameter_covariance=:none, free errors and covariances are NaN, fixed errors remain zero, and diagnostics explain the omission rather than claiming a curvature failure. Conditioning uses free coordinates; a fixed parameter's zero variance is not a singularity. The covariance is a local approximation; profiles and contours remain separate refit operations.

7. Read the result without changing it

report_text, diagnose, and diagnostic_dashboard consume result fields. profile and contour reuse the stored normalized problem, fix one or two parameters, and refit the remaining nuisance parameters. Plotting is an optional CairoMakie package extension consuming the same result objects. None of these paths reruns or alters the original fit unless the API explicitly describes a refit.

Automatic scan ranges intersect the local-covariance range with declared bounds; explicit grids are preserved, and missing threshold crossings remain missing rather than being replaced by a bound, in Gaussian and general-likelihood profiles, contours, and matrices alike.

The Python renderer uses native Matplotlib objects, not Makie. Both renderers share residual/ratio preparation; observation pulls use stored whitened residuals without auxiliary parameter terms, completed profile/contour snapshots render without further scans, and Matplotlib's constrained layout reserves outside legends and reports without a separate layout engine or resize callbacks.

Numerical Invariants

These invariants are the review contract for changes to the numerical core.

Source Map

The core is split by responsibility, not by feature-specific vertical stacks.

This map is also a review rule: a plotting feature should not add a second statistical calculation, and a new optimizer should not own covariance semantics.

Where A New Feature Belongs

Before adding a type or abstraction, ask whether an existing problem can already express the required statistics.

Do not add a parallel result type, cache, or solver path merely to support a new presentation.

Verification Map

Architecture changes need evidence at the layer they affect:

The core gate is julia --project=. test/core_runtests.jl; the complete package gate is julia --project=. -e "using Pkg; Pkg.test()", which adds the optional CairoMakie test dependency. Benchmark methodology is in Performance Checks.

For a page-level output check, run julia --project=docs test/docs_output_snapshots.jl gallery/resonance_decay.md (more page paths are accepted; without arguments it executes every documented workflow). It compares the displayed output with both the page's code cells and the example generator; only solver iteration counts are normalized across solver/platform versions.

Ecosystem Integrations

Packages and Interfaces holds executable examples of distribution fitting, named model construction, interchangeable solvers and posterior inference. The Python interface keeps NumPy callbacks and Matplotlib plots and does not wrap Julia constructor or solver objects.

Analytic references cover normalization, gradients, Hessians, mixture boundaries and nuisance refits. A named extended model agrees between Optim and NativeMinuit; the opt-in probe benchmarks/ecosystem.jl compares both paths at matched accuracy.

Real-Data Reference

The LHCb three-hadron B-decay data supply the mass-spectrum example: 3,420,295 MagnetUp candidates, 9,717 after the notebook's PID cuts, 7,368 in the fit window. examples/data/lhcb_mass/prepare.jl verifies the ROOT checksum and rebuilds every bin; the CC0 collision data, DOI, cuts and model limits are documented on the page, and the full ROOT download stays outside Git and the ordinary docs build.

benchmarks/lhcb_reference.py checks the two peak shapes independently using SciPy CDFs and C++ Minuit2; executed documentation checks the minima, local errors and signal-yield profile against those references. This is a conditional mass-spectrum fit, not a reproduction of an LHCb paper; keep data preparation, warm fitting and uncertainty analysis separate in timing, and do not advertise 3.4 million scanned candidates as 3.4 million fitted events.

The binned distribution adapter processes each mixture component across all bins rather than dispatching per bin, retaining log-space probability sums and derivatives through truncation and zero weights; regressions compare values, gradients and Hessians with the scalar formulation and bound allocations for 10,000 bins. Unbinned heterogeneous mixtures reduce native component batches in blocks of at most 4,096 events, bounding event-sized gradient/Hessian scratch without binning, sampling, or rebuilding the model per block; total work still scales with the event count, and regressions cover the final partial block and zero weights.

Likelihood quadrature estimates its error from the maximum norm of the value and all nested AD coefficients, so a locally constant density cannot hide an oscillatory derivative. Moving finite integration bounds are mapped to [0, 1] without discarding derivatives. CDFs violating range/order by roundoff are reintegrated from the PDF in integration=:auto; strict :cdf and larger violations still fail, and no probabilities are clipped. A DistributionsHEP tail case reproduces the original CDF-roundoff failure.

The example uses upper-only truncation of an exponential's natural support: Distributions 0.25.131's redundant lower bound at zero produces undefined AD derivatives in its log normalizer, while the equivalent upper-only constructor keeps valid automatic derivatives. The shared scalar-solver boundary rejects non-finite initial gradients and, when required, Hessians; the error names the derivative failure and the finite-difference option and never silently changes differentiation policy. No upstream types are patched.

Shared Contract

1. Statistical meaning stays in ScientificFitting. Solvers receive the same validated objective/residuals, parameter mapping, and constraints. Likelihood costs use $-2\log L$; the Minuit adapter must use the matching error scale (errordef=1). Probability densities, discrete masses, and event intensities remain distinct. Dependent observations require a joint likelihood, not a product of marginal probabilities; discrete data do not imply discrete fit parameters are supported.
2. Capabilities are explicit. Check bounds, nonlinear constraints, required derivatives, and residual versus scalar objectives before solving. Reject unsupported requests instead of dropping constraints or adding hidden penalties. Retain explicit solver/derivative options in multistart and all nuisance-parameter refits. Document backend-specific tolerance and budget meanings rather than pretending iterations and function calls are identical.
3. Results preserve evidence. Normalize parameters, objective values, and convergence status while retaining native failure details. Keep minimizer selection independent of local covariance, MINOS intervals, and contour methods. Record missing crossings, parameter boundaries, and failed scans; do not substitute symmetric errors for an invalid asymmetric interval.
4. Metadata is not statistics. Constructor metadata used for optimizer step sizes is not a measurement uncertainty or Gaussian prior. Parameter terms remain explicit. Reject conflicting names/bounds and avoid silently mutating the user's constructor during optimization or profiling.
5. Reuse expensive work at the correct scope. Prepare a parameter-dependent distribution and its normalization once per distinct model at each parameter point, not once per observation. Preserve dual-number types, batched evaluation, stable log densities, and analytic bin integrals/CDF differences where available. Fall back to controlled quadrature or an explicit finite derivative mode where appropriate; never cache a normalization across changed parameters or claim smoothness for an arbitrary interpolated density.

Solver And Posterior Checks

The solver contract is exercised with LsqFit, Optim through Optimization.jl, bounded NLopt Nelder-Mead, and NativeMinuit, across residual, scalar-gradient, nonlinear-constrained, derivative-free and profile-refit paths; C++ Minuit2 provides an independent likelihood-fit reference.

Turing is a complementary posterior-inference workflow, not a MIGRAD substitute. Its external-likelihood interface reuses -problem.objective(p)/2 for an explicitly normalized data likelihood; the executed example specifies priors and support in Turing, without copying SF parameter terms or constraints. A Gamma-Poisson reference checks data/prior separation, values, gradients, Hessians, and four NUTS chains against the analytic posterior and Turing's chain diagnostics. The core has no Turing dependency or sampler wrapper; posterior credible intervals and profile confidence intervals keep their different meanings.

Each required adapter needs analytic or independent parameter/objective/error references, a constraint/failure case, profile-refit parity, and executable documentation. A nested signal-plus-background model must work through BuildConstructors with two solvers without rewriting its likelihood. Measure adapter overhead against direct upstream calls at the same objective and accuracy, separating startup, warm fitting, and uncertainty analysis, by extending existing targeted tests rather than duplicating a benchmark framework.

The core supports Julia 1.10. The optional NativeMinuit 0.7.2 requires Julia 1.11 and declares LGPL-2.1-or-later; the Turing 0.48 reference is executed on the docs environment (Julia 1.13.0). Plot tests cover CairoMakie 0.13 and 0.15.14+, the newer line permitting coexistence with Turing's chain-plotting dependencies.

Performance Checks

The canonical benchmark run is

julia --project=benchmarks benchmarks/runbenchmarks.jl --seconds=1

with --save and --compare for TOML baselines; a comparison fails on missing benchmark cases or mismatched machine metadata. Two standalone probes back the scaling numbers quoted in Validation: benchmarks/whitening_scaling.jl (matrix-free whitening against dense covariance at growing $n$) and benchmarks/ttfx_probe.jl (first-fit and per-model compilation cost); the startup gate below executes benchmarks/startup_probe.jl. Two CI gates back the claims on this page: test/startup_probe_gate.jl verifies fresh-process loading — using ScientificFitting without Makie — together with nested automatic x derivatives, and test/performance_budget_gate.jl checks broad steady-state budgets for representative warmed hot paths, catching regressions such as losing the fast path or recomputing static covariance work inside the objective. The LsqFit least-squares fast path is selected automatically for unbounded static chi-square fits ([-Inf, Inf] counts as no-op bounds), and its weighted Jacobian is reused during FitResult construction.