API Reference

CostFunction{F,T} <: AbstractCostFunction

Wraps a user FCN f::F with an error definition up::T and an internal call counter. Closure-specialized via parametric F (ROADMAP §2.3

• Risk #4) so the call site cf(x) devirtualizes through the concrete

type rather than hitting Julia's ::Function vtable.

Mirrors MnFcn and MnUserFcn from C++ Minuit2, collapsed since Julia doesn't need the C++ inheritance hierarchy (the call counter is in the same struct as the user function; multiple dispatch handles gradient-vs-no-gradient via a separate CostFunctionWithGradient type).

For an unbounded fit the call to f(x) passes the parameter vector unchanged; for a bounded fit the internal→external sin/√ transform is applied on the same call boundary (see transform.jl).

Fields

• f::F — the user function. Must accept an AbstractVector{Float64} (or compatible) and return a Float64-convertible cost value.
• up::T — error definition (1.0 for χ² fits, 0.5 for negative log-likelihood fits; default 1.0). Parametric T keeps the door open for ForwardDiff.Dual{...,Float64} users in Phase 2.1 without re-shuffling the type.
• nfcn::Base.RefValue{Int} — call counter; mutates on each invocation via the call-operator overload. Ref is the idiomatic Julia way to hold mutable state inside an otherwise-immutable struct.
• n_nonfinite::Base.RefValue{Int} — count of calls whose return value was non-finite (NaN/±Inf). Mirrors iminuit's FCN::check_value NaN detection (iminuit src/fcn.cpp), except iminuit warns per occurrence (via MnPrint, suppressed at default print level) while NativeMinuit aggregates here and the MIGRAD drivers warn ONCE at the end of a run that did not end valid (handoff F7/P6; valid fits stay silent). Read via nonfinite_calls.

Performance notes

• This struct is not isbits (because of the Ref). One heap allocation per CostFunction instance. That's fine — we construct one per migrad call.
• The f(x)::Float64 annotation enforces the return type contract at the call boundary; if the user's FCN returns a non-Float64, you'll see a clear runtime error at first call, not silent type instability in the MIGRAD inner loop.

Aliasing contract

The argument vector x passed to your FCN is a borrowed reference to NativeMinuit's internal workspace. The same Vector{Float64} instance is reused across line-search and gradient-calculation calls within a single MIGRAD iteration. Your FCN MUST NOT:

• Retain the vector (e.g. push into a captured array). It will be overwritten by the next call.
• Mutate the vector. Treat it as read-only.

If you need the values later, copy them: xcopy = copy(x).

Examples

julia> cf = CostFunction(x -> sum(abs2, x), 1.0);

julia> cf([1.0, 2.0, 3.0])
14.0

julia> ncalls(cf)
1

julia> reset_ncalls!(cf); ncalls(cf)
0

CostFunctionWithGradient{F,G,T}

Like CostFunction but also carries an analytical gradient function g(x)::Vector{Float64}. Used by migrad's gradient=... keyword to substitute AD-provided or hand-written gradients for the central-difference numerical gradient.

Fields

• f::F — user FCN, f(x::AbstractVector{Float64}) -> Real.
• g::G — gradient function, g(x::AbstractVector{Float64}) -> Vector{Float64}. Must return length-n vector matching length(x).
• up::T — error definition (1.0 for χ², 0.5 for NLL).
• nfcn::Base.RefValue{Int} — FCN call counter.
• ngrad::Base.RefValue{Int} — gradient call counter (separate from nfcn since each gradient call is "one shot" vs the 2n+ FCN calls numerical_gradient! makes per iteration).
• check_gradient::Bool — when true (the default), seed_state validates this gradient against a numerical 2-point estimate at the seed point and warns on disagreement beyond tolerance (the CheckGradient discrepancy check; see [_check_user_gradient]). Mirrors C++ FCNGradientBase::CheckGradient() (default true) gating MnSeedGenerator.cxx:124-144. Pass false to skip the seed-time check (the cross-search probe wrappers do this — the user gradient is already validated at the top-level fit, so re-checking each MINOS/contour re-seed is redundant).

Examples

using ForwardDiff
f = x -> sum(abs2, x .- [1.0, 2.0])
cf = CostFunctionWithGradient(f, x -> ForwardDiff.gradient(f, x))
m = migrad(cf, [0.0, 0.0], [0.1, 0.1])

CostFunctionAD(f, up=1.0; chunk_size=nothing) -> CostFunctionWithGradient

Wrap a user FCN f(x::AbstractVector{<:Real})::Real into a CostFunctionWithGradient whose gradient is computed via ForwardDiff.gradient(f, x).

Requires using ForwardDiff alongside using NativeMinuit (Julia 1.9+ package extension). Without ForwardDiff loaded, calling this throws an informative error.

Arguments

• f — user FCN. Must be generic over element type (i.e. avoid function f(x::Vector{Float64}) ...; write function f(x) ...) so ForwardDiff can promote inputs to ForwardDiff.Dual.
• up::Real=1.0 — error definition (1 for χ², 0.5 for NLL).
• chunk_size::Union{Int,Nothing}=nothing — ForwardDiff chunk size. nothing lets ForwardDiff pick automatically. For n > 12 a smaller chunk (e.g., 4-6) can reduce memory pressure.
• check_gradient::Bool=true — validate the AD gradient against a numerical 2-point estimate at the seed and warn on disagreement (the C++ Minuit2 CheckGradient discrepancy check). Default true matches C++; pass false to skip the one-time seed check.

When to use

• FCN ≥ ~500 ns/call AND n ≤ ~30: AD typically wins over numerical
• HEP amplitude fits with complex intermediates: ForwardDiff works through Complex{Dual}; no special user code needed
• AD-incompatible FCN (mutates Float64 buffers, hard-codes types, calls C libraries): use plain CostFunction instead

Examples

using NativeMinuit, ForwardDiff
chi2(par) = sum((data .- model.(x_data, Ref(par))).^2 ./ errors.^2)
cf = CostFunctionAD(chi2, 1.0)
fmin = migrad(cf, x0, errs)

Beyond C++ Minuit2

C++ Minuit2's MnFcn::operator()(const vector<double>&) is a virtual function — virtual functions cannot be C++ templates, so the input type is locked to double. Generic AD (ForwardDiff Dual, Enzyme, etc.) cannot promote through this signature. Julia has no such barrier: user FCNs are generic on element type by default, so swapping Vector{Float64} for Vector{Dual{...}} requires no user code change. This factory is the visible end of that capability.

migrad(cf, x0, errs;
       strategy=Strategy(0),
       tol=0.1,
       maxfcn=200 + 100·n + 5·n²,
       prec=MachinePrecision()) -> FunctionMinimum

Native-Julia MIGRAD minimization of the user FCN cf starting from x0 with per-parameter step sizes errs. Returns a FunctionMinimum.

This is the low-level unconstrained core — bounds and fixed parameters are handled by the Minuit front end, and analytical/AD gradients via the cost-function type. All Strategy levels 0/1/2 are supported (Strategy ≥ 1 runs the inner-HESSE refinement). Default maxfcn matches C++ MnApplication.cxx:43.

• Convergence: stop when edm ≤ tol · 0.002 (VariableMetricBuilder.cxx:66), or when nfcn ≥ maxfcn.

Arguments

• cf::CostFunction — user FCN wrapper (auto-counts calls).
• x0::AbstractVector{<:Real} — initial parameter values.
• errs::AbstractVector{<:Real} — initial step sizes (≥ 0; algorithm uses |werr| defensively).

Keyword arguments

• strategy::Strategy=Strategy(0) — Strategy level; 0/1/2 all supported.
• tol::Real=0.1 — convergence tolerance on EDM (after *0.002 factor).
• maxfcn::Integer=200+100·n+5·n² — call-count limit.
• prec::MachinePrecision=MachinePrecision() — floating-point precision.

Returns

FunctionMinimum with:

• is_valid=true if MIGRAD converged within tol AND nfcn.
• reached_call_limit=true if nfcn ≥ maxfcn before convergence.
• above_max_edm=true if final EDM > 10·(tol·0.002).
• made_pos_def=true if MnPosDef perturbed the matrix at any point.
• nonfinite_fval=true (⇒ is_valid=false) if the final fval is NaN/±Inf; n_nonfinite_calls counts FCN evaluations that returned a non-finite value during the run. If any occurred AND the fit did not end valid, a single @warn is emitted at the end (suppress with warn_nonfinite=false; valid fits never warn). iminuit parity: NaN FCN values pass through to the minimizer unchanged (iminuit's throw_nan=False default); non-finite trials lose every IEEE < comparison against a finite incumbent, so they are rejected exactly as in C++ Minuit2.

Thread safety

All internal scratch buffers (ping-pong state buffers, line-search scratch, etc.) are stack-local to each migrad call, so multiple threads running migrad concurrently on different CostFunction objects are safe.

CostFunction carries a mutable Base.RefValue{Int} call counter, so sharing one CostFunction across threads causes a benign race on the counter (no memory corruption, but nfcn will be undercounted). Best practice: one CostFunction per thread (cheap — the wrapped function is reused).

migrad(cf::CostFunction, seed::MinimumState;
       strategy=Strategy(0), tol=0.1, maxfcn=..., prec=MachinePrecision())
    -> FunctionMinimum

Low-level MIGRAD entry point that takes a pre-built MinimumState seed. Skips the ~1 + 4·n FCN-call seed_state bootstrap. Intended for callers that already have a warm gradient + invhessian (e.g., the parabolic-fit probe chain in `functioncrossmulti— each probe reuses the previous probe's converged state via [warmrestart_state`](@ref)).

The seed must satisfy is_valid(seed). Caller is responsible for having evaluated the FCN at seed.parameters.x (so the call counter agrees with seed.nfcn); warm_restart_state handles this automatically.

Mirrors the C++ MnApplication::operator()(unsigned int maxfcn, double tolerance) overload that takes an already-constructed MinimumSeed/MinimumState (vs. the user-x0/errs overload which calls MnSeedGenerator).

migrad(f, x0, errs; up=1.0, strategy=Strategy(0), tol=0.1, maxfcn=..., prec=...)

Convenience overload that wraps a bare callable f into a CostFunction before dispatching to the main migrad. The parametric CostFunction{typeof(f)} ensures closure specialization at the call site (parallel-review #2 F4 + ROADMAP §3.4 Criterion 4).

up=1.0 for χ² fits, up=0.5 for negative log-likelihood.

migrad(cf::CostFunctionWithGradient, seed::MinimumState;
       strategy=Strategy(0), tol=0.1, maxfcn=..., prec=..., scratch=nothing)
    -> FunctionMinimum

Phase B's seed-state entry point + Phase F's analytical-gradient dispatch. Same semantics as migrad(::CostFunction, ::MinimumState): skip seed_state's bootstrap, feed the pre-built seed directly into _migrad_loop. Used by warm_restart_state inside MnFunctionCross probes when the user FCN carries an analytical gradient.

migrad(cf::CostFunction, params::Parameters;
       strategy=Strategy(0), tol=0.1, maxfcn=nothing,
       prec=MachinePrecision()) -> BoundedFunctionMinimum

Bound- and fixed-parameter-aware MIGRAD. User FCN receives external parameter values; internal MIGRAD sees only the free parameters in internal (unbounded) coordinates.

Phase 1 first cut. Mirrors the C++ Minuit2 user flow but with the Parameters/Transformation wired in.

Arguments

• cf::CostFunction — wraps the user FCN on EXTERNAL parameter values.
• params::Parameters — parameter metadata: values, errors, bounds, fixed flags, names.

Keyword arguments

• strategy::Strategy=Strategy(0) — Strategy level; 0/1/2 all supported (the Minuit front end defaults to Strategy(1)).
• tol::Real=0.1 — convergence tolerance.
• maxfcn::Union{Integer,Nothing}=nothing — defaults to the standard 200 + 100·n + 5·n² where n is the number of free parameters.
• prec::MachinePrecision.

Returns

BoundedFunctionMinimum with .ext_values and .ext_errors in EXTERNAL coordinates.

migrad(cf::CostFunctionWithGradient, params::Parameters; ...) ->
    BoundedFunctionMinimum

AD-aware bounded MIGRAD. The user FCN + gradient operate on EXTERNAL parameters; this overload threads the int↔ext chain rule into the gradient before handing both to the unbounded migrad(cf::CFwG, ...) path. The 5-10× FCN-call savings of analytical gradients carry through to bounded fits.

Phase F: previously this stored a plain CostFunction view of the internal FCN in BoundedFunctionMinimum.internal_cf, which dropped the analytical gradient on follow-up MINOS / contour calls — codex review caught the regression. The internal CF is now stored as the FULL CostFunctionWithGradient so the AD path survives MIGRAD → MINOS / contour transitions (the downstream function_cross[_multi] / minos / contour_exact signatures already accept AbstractCostFunction).

migrad(m::Minuit; ncall=nothing, resume=true, precision=nothing,
                   strategy=m.strategy, tol=m.tol,
                   iterate=5, use_simplex=false) -> Minuit

IMinuit.jl-compatible alias for migrad!. Mutates m.fmin and returns m. The ncall / resume / precision kwargs are accepted for IMinuit.jl interface parity:

• ncall::Union{Integer,Nothing} ≡ maxfcn cap (default uses NativeMinuit's 200 + 100·n + 5·n² formula).
• resume::Bool=true — if false, reset m.fmin and m.minos_errors before running (matches iminuit's resume argument).
• precision::Union{Real,Nothing} — override the MachinePrecision eps value (rarely used).

strategy and tol default to whatever the user stored on m (settable via m.strategy = ..., m.tol = ...). The constructor default is Strategy(1) — matching iminuit's Minuit-class default and C++ Minuit2's MnStrategy() — for both numerical and analytical/AD (grad=) FCNs. Override per call with strategy=..., or set m.strategy = ... once before the first migrad.

iterate and use_simplex are threaded through to migrad! unchanged. By default (use_simplex=false) the retry is iminuit's _robust_low_level_fit — re-run at the user's strategy, no Simplex, no strategy bump — so this is drop-in-equivalent to iminuit's m.migrad(). The opt-in use_simplex=true enables NativeMinuit's Simplex multistart (which does bump numerical FCNs to Strategy(2)); that is a documented extension beyond C++ Minuit2 / iminuit — see the migrad! docstring.

FunctionMinimum

The result of a MIGRAD minimization. Phase 0 contract:

• state — final MinimumState (parameters, error matrix, gradient, EDM, NFcn).
• seed — the initial state from MnSeedGenerator, kept for reference (e.g. to query starting parameters or initial EDM).
• up — ErrorDef (1.0 for χ², 0.5 for NLL); mirrors the user's CostFunction.up.

Status flags (mirror C++ FunctionMinimum):

• is_valid — overall convergence success.
• reached_call_limit — nfcn ≥ maxfcn was hit.
• above_max_edm — final EDM is more than 10× the requested tolerance.
• hesse_failed — Hesse refinement (Phase 1) failed.
• made_pos_def — MnPosDef perturbed the error matrix.

P6 non-finite-FCN diagnostics (NativeMinuit additions; C++ FunctionMinimum carries no analogue, but the observable verdict matches iminuit — a NaN incumbent there also surfaces as valid=False):

• nonfinite_fval — the final fval is non-finite (NaN/±Inf). Always implies is_valid == false; this flag records the explicit reason. A non-finite fval can only become the incumbent when the FCN is non-finite at the very first evaluation (seed) or when a -Inf trial legitimately wins the IEEE comparisons — mid-run NaN/+Inf trials can never displace a finite incumbent.
• n_nonfinite_calls — how many FCN evaluations during this MIGRAD run returned a non-finite value (0 for a healthy fit). Non-zero with a finite fval means the minimizer brushed an undefined region but the incumbent stayed finite.

M6 (GAP_AUDIT) — per-iteration history:

• states::Vector{MinimumState} — per-iteration snapshots when MIGRAD was invoked with storage_level >= 1. Empty (MinimumState[]) by default (storage_level = 0) — keeps the zero-alloc gate happy.
• storage_level::Int — 0 (default) for no history, 1 for per-iteration snapshots. Mirrors C++ BasicFunctionMinimum.h:109,165.

BoundedFunctionMinimum

Result of bounded migrad(cf, params). Wraps the internal FunctionMinimum and exposes user-visible external values.

Fields

• internal::FunctionMinimum — the MIGRAD result in internal coords.
• params::Parameters — the parameter metadata (bounds, fixed flags).
• ext_values::Vector{Float64} — final external parameter values (length = total n_pars, including fixed parameters at their initial values).
• ext_errors::Vector{Float64} — final external 1σ errors via Jacobian chain rule (length = n_pars; fixed params have 0).
• ext_covariance::Union{Nothing,Matrix{Float64}} — full n_pars × n_pars external covariance matrix (or nothing if inner MIGRAD did not produce a covariance).
• internal_cf::CostFunction — the int-coord-wrapped FCN used by the internal MIGRAD. Required for follow-up calls (MINOS, contour) that operate on internal and must consume internal coordinates — using the user's external cf there would leak coordinate frames (parallel-review #4 A7/B4 blocking).

MinuitParameter

A single fit parameter — name, value, step size, optional bounds, optional fixed flag. Mirrors C++ MinuitParameter from reference/Minuit2_cpp/inc/Minuit2/MinuitParameter.h.

Fields

• name::String — display name.
• value::Float64 — current external value.
• error::Float64 — initial step size (also called "error" in iminuit/Minuit2 nomenclature).
• lower::Float64 — lower bound (NaN if unbounded below).
• upper::Float64 — upper bound (NaN if unbounded above).
• fixed::Bool — true if parameter is currently fixed (excluded from optimization).

NaN is the explicit "absent bound" sentinel, matching the bound_kind classifier in transform.jl.

Parameters

Collection of MinuitParameters plus internal/external index mappings

• precision context. Replaces the tightly-coupled C++ pair

MnUserParameters + MnUserTransformation.

Fields

• pars::Vector{MinuitParameter} — full parameter list (variable + fixed).
• ext_of_int::Vector{Int} — ext_of_int[i_internal] = external_index. Length = number of free parameters.
• int_of_ext::Vector{Int} — int_of_ext[i_external] = internal_index, or 0 if the external parameter is fixed. Length = total parameters.
• name_to_ext::Dict{String,Int} — name → external index (1-based).
• prec::MachinePrecision — used for ext2int clamping in Sin transform.

The mappings are computed once at construction; they don't change as parameters are fixed/released (would require rebuilding).

@enum BoundKind NoBounds BothBounds UpperOnly LowerOnly

Classifies a parameter by which bounds are active. Used by Transformation (in parameters.jl) to select the appropriate int2ext/ext2int/dint2ext functions per parameter.

bound_kind(lower, upper) -> BoundKind

Classify a parameter's bound configuration. NaN is the "absent bound" sentinel.

int2ext(kind, v, lower, upper) -> Float64

Dispatch to the correct internal-to-external transformation by bound kind. For NoBounds, returns v unchanged.

ext2int(kind, ext, lower, upper, prec=MachinePrecision()) -> Float64

Dispatch to the correct external-to-internal transformation by bound kind.

dint2ext(kind, v, lower, upper) -> Float64

Dispatch to d(ext)/d(int). For NoBounds, returns 1.0.

int2ext_error(kind, val, err, lower, upper) -> Float64

External (asymmetric-averaged) parameter error from the internal error err. err is the errordef-SCALED internal 1σ error, i.e. err = sqrt(cov(i,i)) = sqrt(2·up·V_int[i,i]) — NOT the raw sqrt(V_int[i,i]). The C++ comment is err = sigma Value == std::sqrt(cov(i,i)) and the caller passes std::sqrt(2.*up*Error().InvHessian()(i,i)) (reference/Minuit2_cpp/src/MnUserParameterState.cxx:142). Mirrors C++ MnUserTransformation::Int2extError (reference/Minuit2_cpp/src/MnUserTransformation.cxx:115-141).

For unbounded parameters: returns err unchanged.

For bounded parameters: computes the symmetric average of the two-sided perturbations through int2ext:

ui = int2ext(val)
du1 = int2ext(val + err) - ui
du2 = int2ext(val - err) - ui
return 0.5 · (|du1| + |du2|)

with a special clamp for double-bounded when err > 1 (the sin-transform saturates, so |du1| is replaced by the full range).

Phase 1.x D5 (codex parallel-review #4) — closes the near-bound error mis-reporting gap. The Jacobian-diagonal alternative sqrt(V_ext[i,i]) = D · sqrt(V_int[i,i]) underscores near bounds because D = d(ext)/d(int) shrinks toward zero; the two-sided formula captures the actual nonlinear remapping.

hesse(cf, state, strategy=Strategy(1); prec=MachinePrecision(), maxcalls=0)
    -> MinimumState

Compute the full numerical Hessian at state.parameters.x and return a new MinimumState with the refined error matrix, recomputed EDM, and updated FCN call count.

Mirrors C++ MnHesse::operator()(MnFcn, MinimumState, MnUserTransformation, maxcalls) — reference/Minuit2_cpp/src/MnHesse.cxx:93-330.

Arguments

• cf::CostFunction — the user FCN. Operates on the coordinate frame state.parameters.x reports: for a bounded fit this is the INTERNAL (sin/sqrt-transformed) frame supplied by the migrad(cf, params) / hesse(m::Minuit) wrappers, which is exactly the frame the C++ step clamp (has_limits) targets.
• state::MinimumState — current state. The gradient field provides initial step sizes (gst[i] = state.gradient.gstep[i]) and the algorithm refines g2[i].
• strategy::Strategy — controls hessian_ncycles (cycles per parameter), hessian_step_tolerance, hessian_g2_tolerance, and (Strategy ≥ 1) gradient refinement.
• prec::MachinePrecision — floor for step sizes and pos-def gate.
• maxcalls::Integer — FCN call cap; 0 means use the default 200 + 100n + 5n² per MnApplication.cxx:43.
• has_limits::Union{Nothing,AbstractVector{Bool}} — per-parameter (internal index) bound flags. nothing (default) ⇒ all unbounded ⇒ no step clamp ⇒ byte-identical to an unbounded HESSE. When a parameter is flagged, its diagonal probe step d is clamped at 0.5 in internal coordinates (C++ MnHesse.cxx:160-167, 194-195).

Return statuses

The returned MinimumState's error.status may be:

• MnHesseValid — full success.
• MnMadePosDef — pos-def perturbation was applied.
• MnHesseFailed — sag stayed zero or maxcalls hit; matrix is diagonal 1/g2[i] (or 1 where g2 is too small).
• MnInvertFailed — inversion failed; matrix is the same diagonal.

hesse(f, x0, errors; up=1.0, strategy=Strategy(1),
      prec=MachinePrecision(), maxcalls=0, print_level=0) -> HesseResult

Compute a full numerical Hessian (and the derived covariance + errors) for the FCN f at the point x0, using errors as the initial per-parameter step sizes. No MIGRAD is run — this is the standalone MnHesse, the C++ MnHesse::operator()(FCNBase, par, err, maxcalls) overload (reference/Minuit2_cpp/inc/Minuit2/MnHesse.h:57-74).

The errors are computed at the given point, which need not be a minimum — that is the caller's responsibility. At a non-stationary point the Hessian is still well defined; MnPosDef enforcement may flip the status to MnMadePosDef if the curvature is not positive-definite there (e.g. a saddle).

Arguments

• f — the user FCN; takes an AbstractVector{Float64}, returns a real.
• x0::AbstractVector{<:Real} — the point to evaluate at.
• errors::AbstractVector{<:Real} — initial step sizes (one per parameter), the natural scale of each parameter; used to seed the central-difference Hessian and refined internally.

Keyword arguments

• up::Real=1.0 — error definition (1.0 for χ², 0.5 for NLL).
• strategy::Strategy=Strategy(1) — HESSE refinement level.
• prec::MachinePrecision, maxcalls::Integer, print_level::Integer — forwarded to the internal Hessian pass (maxcalls=0 ⇒ default budget).

Returns

A HesseResult exposing .covariance, .errors, .edm, .status, and .valid.

Example

julia> r = hesse(x -> (x[1]-1)^2 + 4(x[2]-2)^2, [0.0, 0.0], [1.0, 1.0]);

julia> r.covariance        # ≈ [1.0 0.0; 0.0 0.25]  (2·up·inv(H), H=diag(2,8))

hesse(m::Minuit; strategy=Strategy(1), maxcall=0) -> Minuit

IMinuit.jl-compatible: re-run a full numerical HESSE at the current converged minimum to refresh the covariance matrix.

Typical use: a fast Strategy(0) MIGRAD leaves the inverse-Hessian as a DFP approximation, which is usually accurate but can drift in ill-conditioned valleys. hesse(m) recomputes the full 2nd-derivative Hessian numerically (mirrors MnHesse invoked standalone in C++) and updates m.fmin in place with the refined ext_covariance and ext_errors.

Strategy(1) is the iminuit default for HESSE; Strategy(2) is more accurate but slower. The maxcall argument is accepted for IMinuit.jl parity but currently unused (the HESSE implementation has its own budget logic).

Internally:

1. Take m.fmin.internal.state (the converged internal-coord state).
2. Call NativeMinuit.hesse(cf_internal, state, strategy) to refresh state.error.inv_hessian via numerical 2nd derivatives.
3. Re-run the same int→ext Jacobian + Int2extError machinery used in migrad(cf, params) to rebuild ext_covariance + ext_errors.
4. Wrap a fresh BoundedFunctionMinimum and overwrite m.fmin.

Returns m for chaining.

covariance(m::FunctionMinimum) -> Symmetric{Float64,Matrix{Float64}}

The covariance matrix 2·up·V where V = inv(H) is the inverse Hessian and up = ErrorDef. For χ² fits with up=1, this is 2·V.

Returns nothing if the minimum doesn't have a valid covariance (!has_covariance(m)).

The matrix is constructed lazily — each call builds a fresh matrix.

eigenvalues(cov::AbstractMatrix{<:Real}) -> Vector{Float64}

Eigenvalues of a covariance matrix, sorted ascending. Mirrors MnEigen::operator() in reference/Minuit2_cpp/src/MnEigen.cxx.

cov may be Symmetric{Float64}, Matrix{Float64}, or any AbstractMatrix{<:Real}. The matrix MUST be symmetric (use Symmetric(cov) to mark it explicitly if needed).

A small or negative eigenvalue (≪ trace / n) is a red flag: the covariance is ill-conditioned, often because two parameters are nearly degenerate. Use global_cc(cov) to identify WHICH parameters are degenerate.

eigenvalues(m::Minuit) -> Union{Nothing,Vector{Float64}}
eigenvalues(bfm::BoundedFunctionMinimum) -> Union{Nothing,Vector{Float64}}

Convenience overloads. Returns nothing if no covariance is available (MIGRAD hasn't run, or it failed to produce one). Otherwise calls eigenvalues on the external covariance matrix.

global_cc(cov::AbstractMatrix{<:Real}) ->
    (cc::Vector{Float64}, valid::Bool)

Global Correlation Coefficient per parameter, mirrors MnGlobalCorrelationCoeff in reference/Minuit2_cpp/src/MnGlobalCorrelationCoeff.cxx.

Returns a pair (cc, valid). cc[i] = sqrt(1 - 1/(C_ii · C⁻¹_ii)) when the denominator is well-conditioned; falls back to 0.0 if the intermediate C_ii · C⁻¹_ii < 1 (numerical near-degeneracy that would otherwise give a NaN). valid = false if the covariance can't be inverted (HESSE didn't converge — rare).

Interpretation:

• cc[i] ≈ 1 → parameter i is strongly determined by the OTHER parameters (statistical degeneracy; consider fixing or re-parametrizing).
• cc[i] ≈ 0 → parameter i is statistically independent of the rest.

Mirrors MnGlobalCorrelationCoeff::MnGlobalCorrelationCoeff.

global_cc(m::Minuit) -> Union{Nothing,Tuple{Vector{Float64},Bool}}
global_cc(bfm::BoundedFunctionMinimum) -> ...

Convenience overloads. Returns nothing if no covariance is available.

CovStatus

Status of a MinimumError covariance matrix. Mirrors the C++ Minuit2 tag types (MnHesseFailed, MnInvertFailed, MnMadePosDef, MnNotPosDef) plus the implicit "valid" case.

Values:

• MnHesseValid — Hesse calculation succeeded; matrix is accurate.
• MnHesseFailed — Hesse failed during refinement.
• MnMadePosDef — Matrix was forced positive-definite via MnPosDef.
• MnInvertFailed — Matrix inversion failed.
• MnNotPosDef — Matrix is not positive-definite.

MnHesseValid is the default for successful MIGRAD-only convergence.

minos(fmin, cf, par_idx; tlr=0.1, maxcalls=1000, sigma=1,
      strategy=Strategy(0), prec=MachinePrecision()) -> MinosError

Compute asymmetric ±σ errors for parameter par_idx. Mirrors MnMinos::Minos(unsigned int, ...) from reference/Minuit2_cpp/src/MnMinos.cxx.

Notes

• Bounded parameters are supported; hitting a bound is a clean termination (recorded in upper_par_limit / lower_par_limit).
• Inner MIGRAD uses Strategy(0) by default. Strategy 1/2 affects the tlr propagation but not HESSE refinement.
• sigma::Real=1 — confidence level in σ-units (P5). Threads up · sigma² into MnFunctionCross's aim (mirrors iminuit's _TemporaryUp). The returned upper / lower then correspond to the k-σ contour on the parameter.

C++ MnMinos algorithm reproduction (X(3872) follow-up)

This function reproduces C++ MnMinos::FindCrossValue (MnMinos.cxx:136-165) including:

1. Linear-correlation pre-shift of OTHER free parameters along the inverse-Hessian direction before the first inner MIGRAD probe (avoids gradient-descent into side basins on non-convex profiles).
2. Outer V → inner prior_cov: the (n-1)×(n-1) minor of state.error.inv_hessian (par_idx row/col removed) is passed to the inner MIGRAD so the DFP starts with the OUTER's correlation structure (mirrors C++ MnMigrad's single-instance covariance reuse).
3. Cold-fallback from warm_state.x on subsequent probes when warm_restart_state fails (negative g2 / edm regression) — the inner DFP rebuilds g2 from scratch but keeps the converged position.
4. ±σ_HESSE placeholder on invalid sides (matches C++ MinosError::Upper() / Lower() at MinosError.h:54).

Known limitations

• Numerical-gradient inner MIGRAD on pathological profiles: with finite-difference gradients (grad= not supplied to Minuit), some strongly-correlated non-convex fits (e.g. X(3872) par[2] / par[3] lower) can converge the inner MIGRAD to the wrong (side) basin even WITH pre-shift + prior_cov, returning the ±σHESSE placeholder. iminuit's numerical-gradient inner MIGRAD finds the same crossings via subtle line-search / step-size differences not yet replicated here. Supplying analytical gradient (grad= AD) closes this gap — see `BenchmarkExamples/X3872dip/benchfull.jl` for a worked example where jmad matches iminuit and jm_num does not. Tracked as a follow-up: needs deeper investigation of Numerical2P step heuristics or Simplex retry within the inner cross-search.

Returns

A MinosError. Use is_valid(e) to check overall success. upper_state / lower_state carry the full parameter vector at the ±σ crossing endpoint (nothing when that side did not converge).

minos(fmin, cf; tlr=0.1, maxcalls=1000, ...) -> Vector{MinosError}

Compute MINOS errors for ALL free parameters. Convenience wrapper that calls the single-parameter overload in turn. Mirrors C++ MnMinos::operator() which iterates over 0..n-1.

minos(m::Minuit, var=nothing; sigma=1, maxcall=0, kwargs...) -> Minuit

IMinuit.jl-compatible alias for minos!. When var is nothing, runs MINOS on all free parameters. var may be an integer index, a String/Symbol name, or a Vector of either.

The sigma kwarg (confidence level in σ-units) is threaded through the MnFunctionCross up · sigma² scaling (P5 — see function_cross for details). At sigma=1 the behavior is C++-MnMinos-identical; at sigma=k the upper/lower errors correspond to the k-σ contour. maxcall (iminuit, singular) caps the FCN calls of each cross-search and tol / toler set the cross-search tolerance; both are forwarded to function_cross via minos!.

MinosError

The asymmetric error result for a single parameter. Mirrors C++ MinosError.

Fields

• par_idx::Int — 1-based parameter index.
• min_par_value::Float64 — the parameter value at the minimum (NOT the function value). Mirrors C++ MinosError::Min() (reference/Minuit2_cpp/inc/Minuit2/MinosError.h:85-86, fMinParValue at line 30-31). Parallel-review #4 B2 — v1 of this field stored fval, which was a semantic mismatch.
• upper::Float64 — upper asymmetric error (x_+σ - x_min).
• lower::Float64 — lower asymmetric error (x_-σ - x_min, ≤ 0).
• upper_valid::Bool, lower_valid::Bool — true if the MINOS analysis completed cleanly on that side. True also when the search saturated against a parameter bound (the corresponding upper_par_limit / lower_par_limit flag is then raised and the published value is the physical bounddistance — `xbound − xmin). This matches iminuit'sm.merrors[name].isvalid` semantics: hitting a bound is a legitimate MINOS termination, not a failure.
• upper_new_min::Bool, lower_new_min::Bool — true if a lower minimum was discovered during the scan (caller should restart MIGRAD from the better point).
• upper_fcn_limit::Bool, lower_fcn_limit::Bool — call budget hit.
• nfcn::Int — total FCN calls across both directions.
• upper_state::Union{Nothing,Vector{Float64}}, lower_state::Union{Nothing,Vector{Float64}} — full parameter snapshot at the ±σ crossing endpoint. nothing when that side did not converge cleanly. Mirrors C++ MinosError::UpperState() / LowerState() (MinosError.h:73-74). Useful for HEP correlated- systematic studies and at-bound diagnostics — see docs/dev/GAP_AUDIT.md M4.

Note on sign convention

upper is positive (one σ to the right), lower is negative (one σ to the left). For a symmetric well-behaved parabolic minimum, upper ≈ -lower ≈ sqrt(2·up·V[i,i]).

minos_upper(m::Minuit, par; kwargs...) -> Float64
minos_lower(m::Minuit, par; kwargs...) -> Float64

Return ONLY the upper (minos_upper, ≥ 0) or lower (minos_lower, ≤ 0) asymmetric MINOS error for parameter par (Integer index or String name). Mirror of C++ MnMinos::Upper / MnMinos::Lower (reference/Minuit2_cpp/inc/Minuit2/MnMinos.h:50-58).

The sign convention matches MinosError.upper / .lower: minos_upper is one σ to the right (positive), minos_lower one σ to the left (negative). The returned value is identical to the corresponding side of a full minos!(m, par) (same function_cross machinery).

Unlike minos!, these are pure queries: they do NOT mutate m (no m.minos_errors update), matching the C++ const accessors. The full asymmetric error is computed internally and the requested side returned — if you need both sides, prefer a single minos!(m, par).

Accepts the same control kwargs as minos!: maxcall, tol / toler, sigma, strategy, print_level.

minos_lower(m::Minuit, par; kwargs...) -> Float64

The lower (≤ 0) side of the asymmetric MINOS error for parameter par (Integer index or String name) — a pure query that does not mutate m. See minos_upper for the full description, sign convention, and the shared control keyword arguments (maxcall, tol/toler, sigma, strategy, print_level).

function_cross(fmin, cf, par_idx, dir; tlr=0.1, maxcalls=1000,
               strategy=Strategy(0), prec=MachinePrecision()) -> MnCross

Find the step multiplier α along parameter par_idx such that the constrained-minimum (other params re-optimized) satisfies f - fmin = up. Used by MINOS (asymmetric errors) and contours.

Arguments

• fmin::FunctionMinimum — the converged MIGRAD result.
• cf::CostFunction — the user FCN (must match the one used for fmin).
• par_idx::Integer — 1-based parameter index to scan along.
• dir::Real — sign of the scan direction (+1.0 for upper error, -1.0 for lower). Combined with the 1-sigma step from state.error.

Keyword arguments

• tlr::Real=0.1 — tolerance. Internal tolerances tlf = tlr·up and tla = tlr mirror C++ MnFunctionCross.cxx:42-44.
• maxcalls::Integer=1000 — call budget across all inner MIGRADs.
• strategy::Strategy=Strategy(0) — passed to inner MIGRADs.
• prec::MachinePrecision.
• sigma::Real=1.0 — confidence level in σ-units. The crossing aim becomes fmin + up · sigma² (mirrors iminuit's minos(cl=) scaling of MnFunctionCross.aim). At sigma=1 the behavior is C++-identical to a single MnFunctionCross call; at sigma=k the returned aopt converges to ≈ k (in the parabolic approximation), so the caller's aopt · σ_1 product is the k-σ error.
• other_param_seed::Union{Nothing,AbstractVector{<:Real}}=nothing — optional length-(n-1) starting vector for the inner MIGRAD on the free parameters (par_idx removed). When provided, used as the COLD- path seed of probe 1; consumed (i.e. cleared to nothing so later probes get warm-state continuation, not the original α=1 seed). Mirrors C++ MnMinos's pre-shift (MnMinos.cxx:143-165) — without it, the inner MIGRAD on strongly-correlated, non-convex profiles can descend into side basins. Caller (minos.jl) supplies the pre- shifted seed; other callers (e.g. contours) leave it nothing.

Returns

MnCross. Check .valid, .new_min, .fcn_limit to interpret.

Known limitations

• par_limit is never raised. Bounded MINOS works correctly (via the internal-coord CF wrap in Minuit.minos! and migrad_bounded.jl), but the boundary-saturation flag isn't surfaced. Equivalent C++ flag: MnCross::CrossParLimit().
• Inner MIGRAD strategy is max(0, strategy.level - 1), matching C++ MnFunctionCross.cxx:106.

MnCross

Result of function_cross. Mirrors C++ MnCross (reference/Minuit2_cpp/inc/Minuit2/MnCross.h).

Fields

• state::MinimumState — the state at the crossing (or current best if invalid).
• aopt::Float64 — the step multiplier at the crossing; NaN if invalid.
• nfcn::Int — cumulative FCN calls made by function_cross.
• valid::Bool — true if a crossing was found within tolerance.
• new_min::Bool — true if a lower minimum was discovered during the scan (Phase 1+ should restart MIGRAD here).
• fcn_limit::Bool — true if the call budget was exhausted.
• par_limit::Bool — true if a parameter bound was hit (Phase 1+ only; always false in first cut).
• ext_state::Union{Nothing,Vector{Float64}} — the inner-bounded-MIGRAD's converged EXTERNAL parameter vector at the crossing (bounded path only; always nothing for unbounded function_cross, and nothing for invalid results). Used by MinosError M4 snapshot fields.

mncontour(m::Minuit, par1, par2; cl=nothing, numpoints=100,
           size=numpoints, kws...) -> Vector{Tuple{Float64,Float64}}

MINOS 2D confidence region boundary — the MnContours profile algorithm (each boundary point re-minimizes all other free parameters via multi-parameter function_cross; not the ellipse approximation). Returns a vector of (x, y) points in physical (external) coordinates — the same frame as m.values, so the points plot directly (for bounded parameters they are mapped back through the sin/√ transform). Mirrors iminuit ≥ 2.0's m.mncontour(par1, par2; cl=...) including its joint- coverage cl semantics:

• cl = nothing (default) → the joint 2-D 68 % confidence region, Δχ² = delta_chisq(0.68, 2) ≈ 2.28 (× m.up).
• 0 < cl < 1 → the joint 2-D region with that probability content.
• cl ≥ 1 → interpreted as nσ: the probability is chisq_cl(cl², 1) (e.g. cl = 1 → 68.27 % → Δχ² ≈ 2.30; cl = 2 → 95.45 % → Δχ² ≈ 6.18).

size is the iminuit-compatible alias for numpoints (it wins if both are passed); sigma is a legacy alias for cl.

Routes through contour_exact with sigma = √(delta_chisq(cl, 2)) — the proper MNCONTOUR algorithm from reference/Minuit2_cpp/src/MnContours.cxx, with the crossing aim scaled to fmin + up·sigma² (mirrors iminuit's temporary-errordef scaling). The faster but approximate ellipse-based contour_ellipse is kept for cases where a quick visual check is enough.

contour_grid(m::Minuit, par1, par2; size=50, bound=2, grid=nothing,
             subtract_min=false) -> ContourGrid

iminuit's Minuit.contour: evaluate the FCN on a 2D grid in the (par1, par2) plane with all other parameters held fixed at their current (best-fit) values — a 2D slice of the FCN, the two-dimensional analogue of profile. No minimization is performed. (Named contour_grid rather than iminuit's contour because the bare name collides with Plots.contour under using NativeMinuit, Plots; this is also what IMinuit.jl exported as contour.)

Arguments

• par1, par2 — Integer index, String, or Symbol name. Must be two distinct free parameters.
• size=50 — number of grid points per axis.
• bound=2 — scan range: a number k means value ± k·σ per axis (σ = current HESSE error), clipped to any parameter limits; or pass ((x_lo, x_hi), (y_lo, y_hi)) explicitly.
• grid=(xs, ys) — explicit grid axes (overrides size/bound).
• subtract_min=false — subtract the grid minimum from the values (fval becomes Δχ²-like).

Returns a ContourGrid; destructures iminuit-style as xs, ys, F = contour_grid(m, "a", "b") with F[i, j] = FCN(xs[i], ys[j]), and plots directly (plot(g) → filled contour; or draw_contour(m, ...)).

ContourGrid

Result of contour_grid — a 2D grid slice of the FCN in the (par_x, par_y) plane with all other parameters fixed at their current (best-fit) values. Mirrors the return of iminuit's Minuit.contour.

Fields

• par_x::Int, par_y::Int — 1-based external parameter indices.
• name_x::String, name_y::String — parameter names (axis labels).
• x::Vector{Float64}, y::Vector{Float64} — the grid axes.
• fval::Matrix{Float64} — FCN values, fval[i, j] = FCN at (x[i], y[j]) (x-major; the plot recipe transposes for Plots' z[j, i] convention).
• value_x::Float64, value_y::Float64 — the central (current/best-fit) parameter values, marked by the plot recipe.
• up::Float64 — the FCN errordef, for user-drawn Δχ² level lines.
• subtracted::Bool — true if the grid minimum was subtracted (subtract_min = true), i.e. fval is already Δχ²-like.

Destructures like iminuit's 3-tuple: xs, ys, F = contour_grid(m, "a", "b").

For confidence regions use mncontour — a grid slice is NOT a confidence region in a multi-parameter fit (see contour_grid).

contour_ellipse(m::Minuit, par_x, par_y; npoints=20, bins=nothing, kwargs...) -> ContoursError

Compute a fast approximate 2D contour — a (possibly asymmetric) error ellipse built from the two axes' MINOS errors + the off-diagonal covariance. par_x / par_y may be Integer or String. The bins=... kwarg is an IMinuit.jl-compatible alias for npoints (takes precedence when both are passed).

Cheap (2 MINOS runs, no boundary search) but assumes the FCN is near-quadratic. For the exact MnContours boundary use mncontour; for an iminuit-style FCN landscape grid use contour_grid. (Named contour before 0.5.0.)

contour_exact(fmin, cf, par_x, par_y; npoints=20, kwargs...) -> ContoursError

C++ Minuit2-equivalent contour algorithm. Replaces the Phase 1 ellipse approximation with the proper boundary search:

1. Compute the 4 MINOS axis crossings (±σx at y=ymin, ±σy at x=xmin).
2. For each subsequent point: find the longest gap in the current boundary, compute its midpoint + perpendicular ray, call function_cross_multi to find the actual boundary along that ray.
3. Insert the new boundary point in the correct position.

Mirrors reference/Minuit2_cpp/src/MnContours.cxx:34-204. Much more accurate than the Phase 1 ellipse approximation for non-quadratic FCNs; slower (each boundary point requires a separate inner MIGRAD chain).

sigma=1 (default) traces the C++-identical boundary FCN = fmin + up — the curve whose axis crossings are the MINOS ±1σ errors and whose joint 2-D coverage is only 39.3 % (James, The Interpretation of Errors, §1.3.3). sigma=k traces FCN = fmin + up·k²; the iminuit-compatible mncontour drives this with k = √(delta_chisq(cl, 2)) to produce joint confidence regions.

ContoursError

Result of contour_ellipse / contour_exact (the boundary-curve contour algorithms). Mirrors C++ ContoursError.

Fields

• par_x::Int, par_y::Int — 1-based parameter indices.
• points::Vector{Tuple{Float64,Float64}} — the contour boundary in (x, y) external coordinates.
• minos_x::MinosError, minos_y::MinosError — the MINOS errors along each axis (used as the basis of the ellipse).
• nfcn::Int — total FCN calls.
• valid::Bool — true if MINOS succeeded on both axes.
• full_points::Vector{Vector{Float64}} — the full n-dim parameter vector at each contour boundary point: the 2 contour coordinates (par_x, par_y) plus the n−2 profiled (re-minimized) values from the inner cross-search, in the same coordinate frame as points. Populated by contour_exact at zero extra cost (the inner re-minimization state is already computed per boundary point — see contour_parameter_sets). Empty for the ellipse approximation contour_ellipse (no inner re-minimization).

profile(m::Minuit, par; bins=100, size=bins, low=0, high=0) ->
    Vector{Tuple{Float64,Float64}}

1D profile of the FCN along par — same as scan but with iminuit's default bins=100 (vs maxsteps=41 for scan). NO inner minimization. Returns (par_value, fval) pairs.

size is the iminuit-compatible alias for bins (review IMPORTANT #3).

Unlike scan, profile is a pure diagnostic — it does NOT move m to the best grid point (no best-value retention, no state change).

mnprofile(m::Minuit, par; bins=30, low=0, high=0) ->
    Vector{Tuple{Float64,Float64}}

1D MINOS profile: at each grid point along par, FIX par at that value and RE-MINIMIZE all other free parameters. Returns (par_value, min_fval) pairs.

This is the "constrained-MIGRAD profile" iminuit exposes as m.mnprofile(par). It's strictly more informative than a bare scan (profile) because it shows the χ² well "as seen by" the nuisance parameters — but each grid point costs an inner-MIGRAD.

Range defaults: if low == high == 0, uses m.values[par] ± 2 · m.errors[par].

Minuit(fcn, x0; names, errors, limits, fixed, up=1.0, prec=...)

iminuit-style wrapper. Constructs the underlying CostFunction and Parameters and exposes mutating MIGRAD / HESSE / MINOS / contour methods plus iminuit-style property access.

Arguments

• fcn — the user function f(x::AbstractVector) -> Real.
• x0::AbstractVector{<:Real} — initial parameter values (external).

Keyword arguments

• names::Vector{<:AbstractString}=["x0", "x1", …, "x(n-1)"] — parameter names (iminuit-style x0, x1, …).
• errors::Vector{<:Real}=fill(0.1, n) — initial step sizes.
• limits::Vector — per-parameter bounds. Each entry may be:
  • nothing for unbounded,
  • (lo, up) for both bounds,
  • (nothing, up) for upper-only,
  • (lo, nothing) for lower-only.
• fixed::Vector{Bool}=fill(false, n).
• up::Real=1.0 — ErrorDef. 1.0 for χ², 0.5 for NLL.
• prec::MachinePrecision.
• threaded_gradient::Union{Bool,Symbol}=false — parallelize the per-coordinate numerical gradient (needs julia -t N). false (default) = serial; true = force threaded, raising ThreadSafetyError if the FCN is not thread-safe; :auto = thread when nthreads()>1 and the FCN probes thread-safe, else warn once and fall back to serial (never throws). The :auto probe runs at most once and is memoized on the fit. No-op for AD (grad=) fits.
• verify_threading::Bool — for threaded_gradient=true, verify thread safety on the first gradient call (default true when forcing threading).

Methods

• migrad!(m; strategy, tol, maxfcn) — run MIGRAD.
• hesse!(m; strategy) — refine the Hessian.
• minos!(m, par_idx_or_name; ...) — single-parameter MINOS.
• minos!(m; ...) — MINOS on all free parameters.
• mncontour(m, par_x, par_y; numpoints) — exact 2D confidence contour; contour_grid (FCN landscape) and contour_ellipse (fast ellipse) for the cheaper variants.

Properties (iminuit-style)

• m.values — external parameter values.
• m.errors — external 1σ errors.
• m.fval, m.edm, m.nfcn, m.valid.
• m.covariance — full external covariance matrix or nothing.
• m.params — the Parameters (name/bounds/fixed structure). Once a fit is cached it reflects the FIT: m.params.pars[i].value/.error equal m.values[i]/m.errors[i] (iminuit parity, issue #38); before any fit (or after reset/a mutation) it holds the constructor-time initial values and steps.
• m.fmin — the underlying BoundedFunctionMinimum (nothing before migrad!).

AbstractFit

Supertype of every NativeMinuit fit-frontend object. Currently the only concrete subtype is Minuit; the abstraction exists so that generic user/ecosystem code can dispatch on f::AbstractFit (the same way IMinuit.jl / iminuit code does — migrad(f::AbstractFit), minos(f::AbstractFit, …)), and so that any future alternative fit frontend can slot in without breaking such code.

IMinuit.jl drop-in note

IMinuit.jl exposes two concrete subtypes — Fit (keyword/scalar-arg fcn(a, b) construction) and ArrayFit (vector fcn(par) construction). That split is a PyCall wrapping artifact: the two need different PyObject construction. NativeMinuit is native Julia and always calls the FCN as f(::AbstractVector) internally (the keyword constructor wraps fcn(a,b) into x -> fcn(x...)), so the two forms have no behavioural difference after construction. We therefore provide Fit and ArrayFit as aliases of Minuit rather than distinct types — a type that doesn't differ in behaviour shouldn't be a distinct type in idiomatic Julia. Code annotating f::Fit / f::ArrayFit / f::AbstractFit or testing f isa Fit keeps working; only code that dispatches Fit and ArrayFit to different methods (rare — the split was never semantic) would need a touch-up.

Fit

IMinuit.jl-compatible alias of Minuit (scalar / keyword-argument construction in IMinuit.jl). In native NativeMinuit it is the same type as Minuit, so existing IMinuit.jl scripts that build Fit(...) work unchanged. See AbstractFit.

ArrayFit

IMinuit.jl-compatible alias of Minuit (vector construction in IMinuit.jl). In native NativeMinuit it is the same type as Minuit, so existing IMinuit.jl scripts that build ArrayFit(...) work unchanged. See AbstractFit.

migrad!(m::Minuit; strategy=m.strategy, tol=m.tol, maxfcn=nothing,
                   iterate=5, use_simplex=false,
                   threaded_gradient=m.threaded_gradient,
                   verify_threading=m.verify_threading,
                   print_level=m.print_level) -> Minuit

Run MIGRAD on m. By default this is drop-in-equivalent to iminuit's m.migrad(): a single MIGRAD followed by iminuit's _robust_low_level_fit retry — if a pass fails to validate (no-improvement / above-max-EDM exit, see C++ VariableMetricBuilder.cxx:278) and the call limit hasn't been reached, re-run MIGRAD from the last converged point at the same strategy (a fresh re-seed, discarding the possibly-degraded DFP inverse-Hessian), up to iterate-1 more times. C++ Minuit2 itself has no retry — this loop is iminuit's addition, reproduced faithfully here. The re-seed lets a stalled fit escape (IAM cold start: S=0 613 → ~383, S=1 330 → ~326 — via iminuit's retry mechanism; the exact basin reached differs from iminuit's on this ill-conditioned problem, see docs/dev/IAMCONVERGENCEGAP.md § Fidelity).

use_simplex=true (opt-in, NOT the default) enables a structured Simplex multistart that is NOT part of C++ Minuit2 or iminuit — a NativeMinuit extension for genuinely multi-minimum landscapes (e.g. the X(3872) J/ψρ + DD̄* multi-dip fit). Each retry pass takes a Nelder-Mead Simplex hop with a geometrically-growing seed step (×1, ×2, ×4, … from the parameter error scale, capped at the physical range) before re-MIGRAD at Strategy(2) for numerical FCNs (the AD path keeps the user strategy). It is this opt-in path's Strategy(2) escalation that walks the IAM xjm WARM start to χ²=322 (PR #10); at the faithful default, xjm converges to iminuit's 325.8 and 322 is reached the C++/iminuit way — by passing strategy=2.

Both paths keep a fixed-point (cycle) detector — re-convergence to an already-visited (x, fval) basin stops the loop (recovering wasted passes on single-basin fits) — and the safety invariant: iterate=N never yields a worse fval than iterate=1 (best_bfm, the lowest-fval pass, is published).

(The opt-in use_simplex=true multistart's multi-scale escape — for genuine multiple local minima in multi-parameter HEP amplitude / phase-shift / LEC fits — does NOT prove the global minimum was found; global optimization is undecidable. The defensible statement: it searches perturbation scales up to re-convergence on a known basin or the physical parameter range.)

iterate=1 disables the retry loop and reproduces single-shot C++-faithful behavior. The safety invariant is guaranteed by construction: iterate=N never yields a worse fval than iterate=1 (the lowest-fval pass is published, tie → the valid one). The number of MIGRAD passes the call executed is recorded in m.n_passes (1 = no retry).

Updates m.fmin. Returns m for chaining. If the constructor was given grad=..., dispatches into the analytical-gradient path on every pass. strategy / tol / print_level / threading flags default to whatever the user stored on m (settable via m.strategy = ..., m.tol = ... or the constructor kwargs).

If a prior m.fmin exists, pass 1 starts from the previous converged point (iminuit-compatible implicit resume). Use reset(m) (or migrad(m; resume=false)) to drop the prior fit and restart from the constructor's initial values. The stored configuration (initial values, user step sizes — _init_params(m)) is NEVER mutated; the carry-forward builds a fresh Parameters only for the duration of the inner MIGRAD call. (The public m.params property reflects the fit once converged — see the Minuit docstring.)

minos!(m::Minuit, par; kwargs...) -> Minuit

Run MINOS for parameter par (integer index or String name). Updates m.minos_errors. Requires m.fmin to be available (call migrad! first). Returns m.

minos!(m::Minuit; kwargs...) -> Minuit

Run MINOS on all free parameters.

hesse!(m::Minuit; kwargs...) -> Minuit

Bang-named alias of hesse: hesse(m) mutates m (it stores the refreshed covariance), so hesse! matches the migrad! / minos! convention. Identical behaviour — use whichever reads better.

HesseResult

Result of a standalone hesse(f, x0, errors) call. Exposes the covariance and errors computed at the supplied point.

Fields

• x::Vector{Float64} — the point the Hessian was evaluated at.
• covariance::Matrix{Float64} — the parameter covariance 2·up·V where V = inv(H) is the inverse Hessian (matches covariance(::FunctionMinimum); for χ² fits with up=1 it is 2·V).
• errors::Vector{Float64} — 1σ errors sqrt(diag(covariance)).
• edm::Float64 — expected distance to minimum at the point.
• nfcn::Int — FCN calls consumed (seed gradient + Hessian passes).
• status::CovStatus — covariance status (see CovStatus).
• valid::Bool — true if the covariance is usable (MnHesseValid / MnMadePosDef).
• state::MinimumState — the full internal state, for advanced consumers (raw inv_hessian, gradient, …).

set_precision(m::Minuit, p::Real) -> Minuit

IMinuit.jl-compatible: override the floating-point precision used by MIGRAD/HESSE/MINOS. The default MachinePrecision() uses 4·eps(Float64) (matching C++ Minuit2's fEpsMac = 4·ε, audit §14); override only when fitting with synthetic-precision FCN models.

is_thread_safe(cf::AbstractCostFunction, x0::AbstractVector;
                errs=fill(0.1, length(x0)), tol=1e-8,
                strategy=Strategy(0), prec=MachinePrecision()) -> Bool

Standalone helper: test whether a user FCN gives the same numerical gradient sequentially vs. threaded at x0. Returns true if safe, false otherwise. Does NOT throw — for the throw version, see the auto-check triggered by migrad(..., threaded_gradient=true, verify_threading=true) (the default for high-level callers).

Use this to probe a new FCN before committing to threaded_gradient=true:

using NativeMinuit
cf = CostFunction(my_chi2)
if Threads.nthreads() > 1 && NativeMinuit.is_thread_safe(cf, x0)
    m = Minuit(my_chi2, x0; threaded_gradient=true)
else
    m = Minuit(my_chi2, x0)
end

Cost: equivalent to 2 gradient evaluations (≈ 4·n·grad_ncycles FCN calls).

ThreadSafetyError(message)

Raised by _migrad_loop when threaded_gradient=true is used with verify_threading=true and the user FCN's threaded gradient does not match its sequential counterpart at the seed point (Phase H safety check). Indicates the FCN has hidden mutable state that races under parallel evaluation.

simplex(cf::CostFunction, x0, errs;
        maxfcn=nothing, minedm=nothing, prec=MachinePrecision()) -> FunctionMinimum

Nelder-Mead simplex minimizer. Gradient-free; calls cf(x) only. Mirrors SimplexBuilder::Minimum in reference/Minuit2_cpp/src/SimplexBuilder.cxx.

Arguments

• cf::CostFunction — user FCN wrapper. cf.up is the ErrorDef (1.0 for χ², 0.5 for NLL) used to scale the final per-parameter errors.
• x0::AbstractVector{<:Real} — initial parameter values.
• errs::AbstractVector{<:Real} — initial 1σ step sizes. The simplex initial edge length is 10·Gstep[i] per dim with the seed Gstep[i] = max(gsmin, 0.1·|errs[i]|), i.e. an effective edge ≈ |errs[i]| (matches the C++ default step = 10 · seed.Gradient().Gstep(); SimplexBuilder.cxx:38 + InitialGradientCalculator.cxx:64).

Keyword arguments

• maxfcn::Union{Integer,Nothing}=nothing — FCN call budget; defaults to 200 + 100·n + 5·n² (same as MIGRAD).
• minedm::Union{Real,Nothing}=nothing — convergence tolerance on edm = f(jh) - f(jl). Defaults to 0.1 · up, the C++/iminuit Simplex EDM goal: ModularFunctionMinimizer::Minimize scales effective_toler = toler·Up() with the canonical toler = 0.1 and passes it to ALL builders (ModularFunctionMinimizer.cxx:175). The extra ×0.002 of VariableMetricBuilder.cxx:66 is MIGRAD-only — Simplex does not apply it. Pass a smaller minedm explicitly for a tighter (non-C++-default) stopping rule.
• prec::MachinePrecision — floating-point precision (rarely tuned).
• warn_nonfinite::Bool = true — emit the single end-of-run warning when the FCN returned non-finite values AND the run did not end valid (P6, same policy as migrad). Inner probes (the use_simplex multistart inside migrad!) pass false.

Returns

FunctionMinimum with state.parameters.x = best vertex, errors derived from the simplex extent (no inverse Hessian — covariance is nothing). Run a follow-up hesse(cf, state) if you need a covariance.

A non-finite final fval is never valid (P6): the result comes back is_valid = false with the explicit nonfinite_fval = true reason, and n_nonfinite_calls counts the FCN evaluations that returned NaN/±Inf during this run. iminuit parity: m.simplex() on an all-NaN FCN reports fval=nan, valid=False.

Notes

• This is a direct port of the C++ Nelder-Mead recipe (reflect / expand / contract / shrink). The "rho" adaptive step is the L162-176 block — a parabolic fit through (y_h, y_star, y_stst) choosing the next probe location.
• Without a gradient, EDM is not the variable-metric estimator (gᵀV g) but the simpler f_high - f_low over the simplex.

simplex(cf::CostFunction, params::Parameters; ...) -> BoundedFunctionMinimum

Bound-aware Nelder-Mead. User FCN receives external coordinates; the inner simplex runs in internal (transformed) coordinates so bounds and fixed parameters are respected.

Like its unbounded sibling, this overload does not produce a covariance matrix (simplex has no inverse Hessian). m.matrix, eigenvalues(m), global_cc(m) all return nothing after a simplex-only fit — run hesse on top if you need a real cov. External per-parameter errors are derived from the simplex dirin (simplex extent × √(up/edm)) via the int→ext Jacobian.

simplex(m::Minuit; maxfcn=nothing, ncall=maxfcn, minedm=nothing) -> Minuit

Run Nelder-Mead simplex on m. Updates m.fmin and returns m. Useful as a robust fallback when MIGRAD fails (no gradient needed).

ncall is the iminuit-compatible alias for maxfcn (review IMPORTANT #3).

scan(cf::CostFunction, x0, errs, par_idx;
     maxsteps=41, low=0.0, high=0.0) ->
    Vector{Tuple{Float64,Float64}}

1D parameter scan. Evaluates cf(x) at maxsteps equally-spaced points along parameter par_idx, holding other parameters fixed at the values in x0.

If low == high == 0, the scan range defaults to x0[par_idx] ± 2·errs[par_idx].

Returns a vector of (parameter_value, fcn_value) pairs. The first entry is always the central point (x0[par_idx], cf(x0)) — useful for assessing how far the scan extends from the minimum. The remaining maxsteps entries are the equally-spaced scan points.

Mirrors MnParameterScan::operator() in reference/Minuit2_cpp/src/MnParameterScan.cxx.

scan(cf::CostFunction, params::Parameters, par_idx;
     maxsteps=41, low=0.0, high=0.0) -> Vector{Tuple{Float64,Float64}}

Bound-aware 1D scan. Operates in EXTERNAL coordinates (user FCN sees physical values). If the scan range straddles a parameter bound, the range is clipped to the bound. If par_idx is fixed, throws an error.

scan(m::Minuit, par; maxsteps=41, low=0, high=0) ->
    Vector{Tuple{Float64,Float64}}

1D scan along parameter par (Integer index or String name) of the fitted Minuit. Returns (x, fval) pairs; the central point comes first followed by maxsteps equally-spaced probes.

If low == high == 0, defaults to ±2σ around the current parameter value (or the parameter limits if both are set). The scan does NOT minimize over other parameters — for that use mnprofile.

Best-value retention (iminuit / C++ MnParameterScan semantics): the scan runs around the CURRENT values (after a fit, the held parameters sit at m.fmin's converged values — what m.values/m.params now report — not the constructor initials), and as a side effect m is left at the lowest-fval grid point found (the central point is included in the comparison). So m.values / m.fval reflect the best grid point — with the other parameters' current values preserved — and a follow-up migrad!/hesse resumes from there. The covariance is NOT updated (scan computes no Hessian; m.matrix → nothing). The returned point-list is unchanged. To scan WITHOUT moving m, use profile.

fix!(m::Minuit, par::Union{Integer,AbstractString}) -> Minuit

Mark parameter par as fixed (excluded from optimization). Mirrors C++ MnUserParameters::Fix(i) / Fix(name). Drops m.fmin and clears m.minos_errors. Returns m for chaining.

par is either the 1-based external index or the parameter name. Already-fixed parameters are still re-fixed (no-op on the fixed flag, but cache is still invalidated for consistency).

fix!(m, 1)
fix!(m, "alpha")

release!(m::Minuit, par::Union{Integer,AbstractString}) -> Minuit

Clear parameter par's fixed flag (re-include in optimization). Mirrors C++ MnUserParameters::Release(i) / Release(name). Drops m.fmin and clears m.minos_errors. Returns m for chaining.

Pairs with fix! for fix-fit-release-fit profile-likelihood scans:

fix!(m, "alpha"); migrad!(m); release!(m, "alpha"); migrad!(m)

set_value!(m::Minuit, par::Union{Integer,AbstractString}, v::Real) -> Minuit

Set parameter par's initial value to v. Mirrors C++ MnUserParameters::SetValue(i, v) / SetValue(name, v). Drops m.fmin and clears m.minos_errors. Returns m for chaining.

v must be finite (NaN / ±Inf throw ArgumentError) — matches the iminuit Python wrapper's setattr guard. The int↔ext transform clamps to bounds at minimization time, so no value-vs-limits check is done here (matches C++).

set_error!(m::Minuit, par::Union{Integer,AbstractString}, e::Real) -> Minuit

Set parameter par's step size to e. Mirrors C++ MnUserParameters::SetError(i, e) / SetError(name, e). Drops m.fmin and clears m.minos_errors. Returns m for chaining.

e must be finite and non-negative (NaN / ±Inf / negative throw ArgumentError). A non-positive step would degrade the numerical gradient floor and could silently poison the seed.

set_limits!(m::Minuit, par, lo, up) -> Minuit

Set parameter par's bounds to (lo, up). Mirrors C++ MnUserParameters::SetLimits(i, lo, up) / SetLimits(name, lo, up). Drops m.fmin and clears m.minos_errors. Returns m for chaining.

lo and up may each be a Real, nothing, or ±Inf — the latter two are stored as NaN (the "absent bound" sentinel), so passing (nothing, 10.0) makes par upper-bounded only, and (nothing, nothing) is equivalent to remove_limits!. When both are finite Reals, lo < up is required (else ArgumentError).

set_lower_limit!(m::Minuit, par::Union{Integer,AbstractString}, lo::Real) -> Minuit

Constrain par from BELOW ONLY: set its lower bound to lo and clear any upper bound, leaving it half-open [lo, +∞). Mirrors C++ MnUserParameters::SetLowerLimit → MinuitParameter::SetLowerLimit (reference/Minuit2_cpp/inc/Minuit2/MinuitParameter.h:131-137), which sets fLoLimValid=true, fUpLimValid=false. Drops m.fmin and clears m.minos_errors. Returns m for chaining.

NB: clearing the upper bound is the C++ behavior — to instead KEEP an existing upper bound and add a lower one, use the two-sided set_limits! (or m.limits[par] = (lo, hi)). lo must be finite (NaN / ±Inf throw ArgumentError); use remove_limits! to drop a bound.

set_upper_limit!(m::Minuit, par::Union{Integer,AbstractString}, hi::Real) -> Minuit

Constrain par from ABOVE ONLY: set its upper bound to hi and clear any lower bound, leaving it half-open (-∞, hi]. Mirrors C++ MnUserParameters::SetUpperLimit → MinuitParameter::SetUpperLimit (reference/Minuit2_cpp/inc/Minuit2/MinuitParameter.h:123-129), which sets fLoLimValid=false, fUpLimValid=true. Drops m.fmin and clears m.minos_errors. Returns m for chaining.

NB: clearing the lower bound is the C++ behavior — to instead KEEP an existing lower bound and add an upper one, use the two-sided set_limits! (or m.limits[par] = (lo, hi)). hi must be finite (NaN / ±Inf throw ArgumentError); use remove_limits! to drop a bound.

remove_limits!(m::Minuit, par::Union{Integer,AbstractString}) -> Minuit

Clear both bounds on parameter par. Mirrors C++ MnUserParameters::RemoveLimits(i) / RemoveLimits(name). Drops m.fmin and clears m.minos_errors. Returns m for chaining.

After the call, m.params.pars[par].lower and .upper are both NaN.

Data(x, y, err)

Holds x, y, and symmetric err-on-y vectors for χ² fitting. Fields x, y, err, ndata. Mirrors IMinuit.jl's Data struct.

Different Data sets can be concatenated with vcat(d1, d2, ...) and sliced with d[idx]. Asymmetric errors are not supported (use a custom FCN with two separate σ for that).

chisq(dist::Function, data::Data, par; fitrange=()) -> Float64
chisq(dist::Function, data, par;       fitrange=()) -> Float64

χ² cost: dist(x, par) evaluates the model at scalar x with the parameter container par. Returns $\sum_i ((y_i - \mathrm{dist}(x_i, par)) / \sigma_i)^2$.

The second method accepts data as a tuple (x, y, err) or (x, y) (unit errors). fitrange restricts the sum to a subset of data indices (default: all points).

Use chisq as the FCN argument to Minuit/migrad, or wrap via model_fit / @model_fit.

model_fit(model::Function, data::Data, start_values; kws...) -> Minuit
model_fit(model::Function, data::Data, fit::Minuit;    kws...) -> Minuit

Build a Minuit fit of model(x, par) against data using the χ² cost from chisq. start_values may be an AbstractVector of initial values or a previous Minuit fit (whose latest values are reused as the new starting point).

kws... flow through to the Minuit constructor (name, error, limits, fixed, grad, strategy, tol, etc.).

@model_fit model data start_values kws...

Macro form of model_fit. Expands to Minuit(par -> chisq(model, data, par), start_values; kws...).

args(m::Minuit) -> Vector{Float64}

IMinuit.jl-compatible convenience: returns the current parameter values as a Vector{Float64}. Equivalent to m.values.

matrix(m::Minuit; correlation=false, skip_fixed=true) -> Matrix{Float64}

IMinuit.jl-compatible covariance matrix accessor.

• correlation=false (default): returns the external covariance.
• correlation=true: returns the correlation matrix C[i,j] = V[i,j] / √(V[i,i]·V[j,j]).
• skip_fixed=true (default): returns the nfree × nfree submatrix (the free_covariance shape, matching C++ MnUserParameterState).
• skip_fixed=false: returns the full ntotal × ntotal matrix with zero rows + cols for fixed parameters.

Returns nothing if MIGRAD hasn't been called or the covariance is unavailable.

chi2(y, yerror, ymodel) -> Float64

Pearson χ² for symmetric Gaussian errors. Mirrors iminuit.cost.chi2(y, yerror, ymodel). Computes $\sum_i ((y_i - y_{model,i}) / \sigma_i)^2$ skipping bins with yerror[i] ≤ 0.

poisson_chi2(n, mu) -> Float64

Likelihood-ratio χ² for Poisson-distributed counts n with predicted means mu. Mirrors iminuit.cost.poisson_chi2(n, mu): $2 \sum_i \left[ \mu_i - n_i + n_i \log(n_i / \mu_i) \right]$ with the n_i log(n_i / μ_i) term defined as 0 when n_i = 0 (limit of x log x at 0).

multinominal_chi2(n, mu) -> Float64

Likelihood-ratio χ² for multinomial-distributed counts (only relative proportions matter — constant offsets cancel). Mirrors iminuit.cost.multinominal_chi2(n, mu): $2 \sum_i n_i \log(n_i / \mu_i)$ with the 0 log 0 = 0 convention.

func_argnames(f::Function) -> Vector{Symbol}

Argument names of f (skipping the implicit f::Function slot). Useful for auto-deriving parameter names from an FCN signature (func_argnames(par_a, par_b) -> [:par_a, :par_b]).

Strategy(level::Integer)

Mirror of MnStrategy from reference/Minuit2_cpp/src/MnStrategy.cxx:33–70.

Controls the cost/precision trade-off in gradient and Hessian calculations via seven tunable knobs. Three preset levels:

• Strategy(0) — Low: fastest, lowest accuracy. The default for the low-level migrad(cf, …) / seed / function_cross entry points (pinned to the C++ oracle reference data).
• Strategy(1) — Medium: matches the C++ Minuit2 MnStrategy() default and the iminuit Minuit class default. The high-level Minuit(fcn, x0) constructor default (so a bare migrad!(m) is drop-in-equivalent to iminuit's m.migrad()).
• Strategy(2) — High: slowest, highest accuracy.

The level field exposes the integer level (0/1/2) — VariableMetricBuilder.cxx branches on Strategy() >= 1 to invoke the inner MnHesse refinement path. All three levels are supported (hesse.jl shipped in Phase 1); the inner-HESSE re-seed at Strategy ≥ 1 is what lets the DFP loop reach deeper minima on stiff fits — see docs/dev/IAM_CONVERGENCE_GAP.md and docs/dev/DESIGN.md DR-008.

Like MachinePrecision, this is an immutable, isbits, type-stable struct — zero allocation when passed through the call chain.

Examples

julia> s = Strategy(0);

julia> s.level, s.grad_ncycles, s.grad_step_tolerance
(0, 2, 0.5)

julia> Strategy(3)
ERROR: ArgumentError: ...

MachinePrecision(eps_value = 4 * eps(Float64))

Mirror of MnMachinePrecision from reference/Minuit2_cpp/inc/Minuit2/MnMachinePrecision.h.

• eps — relative floating-point precision (C++ fEpsMac). Defaults to 4·eps(Float64) ≈ 8.88e-16, matching C++ MnMachinePrecision.cxx:26 (fEpsMac = 4. * numeric_limits<double>::epsilon()). The ×4 is the C++ factor that absorbs the 2× between numeric_limits::epsilon and the DLAMCH-style epsilon Minuit2 was tuned against (see the C++ source note).
• eps2 — 2·√eps (C++ fEpsMa2), the tolerance multiplier used by the numerical-gradient step-size algorithm. With the default eps this is 2·√(4·eps(Float64)) = 4·√eps(Float64) ≈ 5.96e-8.

A user can override the default to declare reduced precision when the FCN value is itself computed at less than IEEE-754 nominal accuracy (e.g. stochastic simulation FCNs).

This is an immutable, isbits, type-stable struct: stack-allocated when used as a local, no GC pressure in hot loops.

Examples

julia> p = MachinePrecision();

julia> p.eps == 4 * eps(Float64)
true

julia> p.eps2 ≈ 2 * sqrt(4 * eps(Float64))
true

julia> p_noisy = MachinePrecision(1e-12);  # FCN noisy to 1e-12

julia> p_noisy.eps2 ≈ 2 * sqrt(1e-12)
true

AbstractCost

Supertype of the Julia-native cost-function family — LeastSquares, UnbinnedNLL, BinnedNLL, ExtendedUnbinnedNLL, ExtendedBinnedNLL, and their CostSum composition.

A cost is a callable mapping a parameter vector to a scalar objective: (c::AbstractCost)(par) returns a Float64 on ordinary numeric calls, and a ForwardDiff Dual under automatic differentiation (every cost in this family is AD-generic, given an AD-able user model). Two things make a cost more than a bare closure:

• errordef(c) — a trait (dispatched on the concrete type) giving the Minuit ErrorDef (1.0 for χ², 0.5 for a negative-log-likelihood). Minuit(c, x0) reads it automatically, so you never pass up=.
• a + b composes costs into a CostSum for simultaneous fits, taking the union of their (named) parameters.

Distinct from AbstractCostFunction (fcn.jl), which is the internal FCN wrapper (CostFunction / CostFunctionWithGradient) that the MIGRAD/HESSE/MINOS machinery consumes. A cost becomes an FCN when you hand it to Minuit: it is wrapped in a CostFunction carrying up = errordef(c).

LeastSquares(x, y, yerror, model; mask=nothing, name=nothing)
LeastSquares(data::Data, model;   mask=nothing, name=nothing)

Least-squares cost $\sum_i ((y_i - \mathrm{model}(x_i, par)) / \sigma_i)^2$ (errordef 1.0). model(x_scalar, par) follows the NativeMinuit / IMinuit.jl convention (the same model you pass to chisq / model_fit).

This is a thin object-style surface over the shared χ² kernel: an unmasked LeastSquares(x,y,ye,model)(par) is bit-identical to chisq(model, Data(x,y,ye), par). Both are just _chisq_core.

mask is a BitVector (true keeps the point) applied with no data copy. name (a vector of names) is only needed to compose this cost with +; Minuit(LeastSquares(...), x0) otherwise names parameters x0, x1, … like the bare-FCN constructor.

Examples

julia> c = LeastSquares([0,1,2], [1,3,5], [0.1,0.1,0.1], (x,p)->p[1]*x+p[2]);

julia> c([2.0, 1.0])         # perfect fit ⇒ χ² = 0
0.0

julia> m = Minuit(c, [1.0, 0.0]); migrad!(m); m.values
2-element ...:
 2.0
 1.0

UnbinnedNLL(x, pdf; log=false, mask=nothing, name=nothing)

Unbinned negative-log-likelihood $-\sum_i \log \mathrm{pdf}(x_i, par)$ (errordef 0.5). pdf(x_scalar, par) must be normalised over the observable range. Pass log=true if your function returns the log-density directly (numerically safer in the tails).

mask (true keeps the sample) and name behave as in LeastSquares.

ExtendedUnbinnedNLL(x, density, integral; log=false, mask=nothing, name=nothing)

Extended unbinned negative-log-likelihood $\mu(par) - \sum_i \log \rho(x_i, par)$ (errordef 0.5), where density is the differential intensity $\rho = \mathrm{d}N/\mathrm{d}x$ (density(x_scalar, par)) and integral(par) is its integral over the observable range — the expected total event count $\mu$. The extended likelihood fits the normalisation as well as the shape.

This is the scalar-convention split of iminuit's scaled_pdf, which returns (integral, density_array) from one call. Pass log=true if density returns its own logarithm.

BinnedNLL(n, xe, cdf; mask=nothing, name=nothing)

Binned negative-log-likelihood for histogram counts n (length nbins) with bin edges xe (length nbins+1). cdf(x_scalar, par) is the normalised cumulative distribution; the per-bin probability is $p_i = \mathrm{cdf}(x_{i+1}) - \mathrm{cdf}(x_i)$, and the fit uses the multinomial likelihood (conditioning on the observed total). errordef 0.5. Use ExtendedBinnedNLL to also fit the normalisation.

The value equals 0.5 * multinominal_chi2(n, μ) (the existing likelihood-ratio kernel), with μ scaled to the kept-bin total.

mask (true keeps the bin) and name behave as in LeastSquares.

ExtendedBinnedNLL(n, xe, scaled_cdf; mask=nothing, name=nothing)

Extended binned negative-log-likelihood: like BinnedNLL but scaled_cdf(x_scalar, par) returns the expected cumulative count (its integral is the expected total), so the per-bin expectation is $\mu_i = \mathrm{scaled\_cdf}(x_{i+1}) - \mathrm{scaled\_cdf}(x_i)$ and the normalisation is fitted. Poisson likelihood; errordef 0.5.

The value equals 0.5 * poisson_chi2(n, μ).

CostSum(costs...)
c1 + c2

Sum of cost functions for a simultaneous fit. The parameter set is the union (by name) of the components, so components sharing a parameter name share that parameter. Works for mixed cost types (LeastSquares + UnbinnedNLL).

Every component must carry parameter names (name= at construction); the union becomes the default parameter names of Minuit(c1+c2, x0). errordef is 1.0: each component is rescaled by 1/errordef onto the common χ² scale, so the combined fit's uncertainties match a single joint likelihood (and iminuit).

errordef(cf::CostFunction)

Return the error definition (up) — 1.0 for χ², 0.5 for NLL.

errordef(c::AbstractCost) -> Float64

Minuit ErrorDef for the cost: 1.0 for χ² (LeastSquares, CostSum), 0.5 for the likelihood costs. Pure trait dispatch on the type.

Minuit(c, x0) folds this into the underlying CostFunction.up, which HESSE (2·up·H⁻¹) and MINOS then use to scale the parameter uncertainties. (errordef is also defined on CostFunction in fcn.jl, where it returns the stored up.)

delta_chisq(cl, ndof) -> Float64

The Δχ² threshold χ²(ndof)-quantile at confidence level cl, i.e. the amount by which χ² rises from its minimum at the edge of the cl confidence region for ndof jointly-estimated parameters.

cl convention (matches iminuit)

• 0 < cl < 1 — interpret cl as a probability (e.g. 0.95).
• cl ≥ 1 — interpret cl as nσ (Gaussian-equivalent): 1→68.27 %, 2→95.45 %, 3→99.73 %.

⚠ Joint vs. single-parameter — READ THIS

ndof is the number of parameters defining the region, NOT the total number of fit parameters. The two questions below give different Δχ²:

A common mistake is to use Δχ²=1 for a 2-D contour: the 68 % joint 2-parameter region is Δχ²=2.30, not 1. The Monte-Carlo sampler get_contours_samples samples all free parameters jointly, so its default threshold uses ndof = n_free.

Examples

delta_chisq(0.6827, 1)  # ≈ 1.00   (1σ, one parameter)
delta_chisq(0.6827, 2)  # ≈ 2.30   (1σ joint, two parameters)
delta_chisq(0.6827, 3)  # ≈ 3.53   (1σ joint, three parameters)
delta_chisq(0.95,   1)  # ≈ 3.84   (95 %, one parameter)
delta_chisq(1,      2)  # ≈ 2.30   (cl=1 ⇒ 1σ ⇒ 0.6827 ⇒ 2.30)
delta_chisq(2,      1)  # ≈ 4.00   (2σ, one parameter)

See also chisq_cl (the inverse), and docs/src/error_analysis.md.

chisq_cl(dchisq, ndof) -> Float64

Inverse of delta_chisq: the probability (confidence level, in (0,1)) that a χ²(ndof) variate is ≤ dchisq. Equivalently the χ²(ndof) CDF at dchisq.

chisq_cl(1.0,  1)  # ≈ 0.6827
chisq_cl(2.30, 2)  # ≈ 0.6827
chisq_cl(3.84, 1)  # ≈ 0.95

Round-trips with delta_chisq when cl is a probability: chisq_cl(delta_chisq(p, k), k) ≈ p for 0 < p < 1.

get_contours_samples(m::Minuit; kwargs...) -> NamedTuple
get_contours_samples(m::Minuit, χsq, paras=nothing, ranges=nothing; kwargs...)

Monte-Carlo true-Δχ² error region for a converged (or even an ill-converged) Minuit fit — the principled alternative to MINOS when the χ² posterior is non-parabolic or m.fmin is unreliable.

The method (mirroring the X(3872) published analysis):

1. Draw nsamples proposals for all free parameters jointly from MvNormal(best, inflate²·Σ) (proposal=:mvnormal, the default) or uniformly over a user box (proposal=:uniform with ranges).
2. Keep a proposal iff its TRUE Δχ² satisfies χ²(x) − χ²_min ≤ up · delta_chisq(cl, ndof) — the exact FCN is re-evaluated at every sample; the Gaussian is ONLY the proposal, never the acceptance test.
3. Report the kept sets and their per-parameter (min, max) extents.

Proposal under-coverage (the critical pitfall)

A MvNormal(best, Σ) proposal severely under-estimates the region when Σ is unreliable (poorly-converged / invalid fit, or MnPosDef-forced covariance) or when the true region is highly nonlinear and extends beyond the local Gaussian. This routine guards against that:

• inflate::Real=1 — sample MvNormal(best, inflate²·Σ) (manual widening).
• adaptive::Bool=true — detect proposal-limited acceptance (the accepted extent filling the proposed extent) and geometrically grow inflate (widen_factor, up to max_widen_rounds), re-sampling until the region is covered. widen_rounds and under_coverage are reported.
• proposal=:uniform + ranges — a covariance-FREE box proposal that does not depend on Σ at all (use when Σ is meaningless).
• On !is_valid / made_pos_def / hesse_failed, a MvNormal run emits a warning recommending these mitigations rather than silently under-estimating.

Keyword arguments

• χsq — optional explicit χsq(x_full)::Real over the full external parameter vector. Defaults to the fit's own FCN (m.fcn).
• paras — parameters (indices or names) to report bounds for; default all free parameters. (Sampling is always over all free parameters.)
• nsamples::Integer=10_000 — proposals per round.
• cl::Real=1 — confidence level (<1 probability, ≥1 nσ; see delta_chisq).
• ndof::Integer — dof for the Δχ² threshold; default n_free (the joint region over all sampled parameters). Keep ndof = n_free unless you know you want a sub-dimensional threshold.
• proposal::Symbol=:mvnormal — :mvnormal or :uniform.
• ranges — (lo,hi) per free parameter; required for :uniform.
• inflate, adaptive, max_widen_rounds, widen_factor — see above.
• clip_threshold::Real=0.9 — a round is "proposal-limited" when the accepted samples fill more than this fraction of the proposed extent in some parameter; triggers widening.
• extent_tol::Real=0.02 — the accepted region is "captured" (stop widening) once its per-parameter (min,max) change between rounds drops below this relative tolerance.
• threaded::Bool=false — evaluate χ² over samples in parallel. Phase-H aware: the FCN thread-safety is verified first (on a throwaway counter so m.nfcn is untouched); a racey FCN falls back to serial with a warning (a racing FCN would corrupt acceptance). The χ² evaluations bypass the fit's call counter, so sampling does not change m.nfcn. Bounded parameters are sampled in external coordinates; proposals outside a parameter's limits are rejected before the FCN is called.
• seed::Integer — RNG seed for reproducible proposals (proposals are drawn serially, so results are reproducible regardless of threaded).
• mahalanobis::Bool=false — also return the per-sample Mahalanobis distance (DIAGNOSTIC only — never used for acceptance).

Returns

A NamedTuple with fields: samples (n_accepted × n_free matrix, one row per kept set), bounds (Vector{(min,max)} per reported parameter), best, names, free_names, delta_chisq_values (per-sample true Δχ²), n_accepted, n_total, acceptance (kept fraction — for a Gaussian posterior this is ≈ cl only when ndof == n_free, a useful sanity check), widen_rounds, inflate_final, delta (the Δχ² threshold in χ² units), up, cl, ndof, proposal, under_coverage (true if widening hit its cap, or stopped on a still-clipping round, while the region was plausibly larger than sampled), and mahalanobis (or nothing).

See contour_df_samples for a DataFrame of samples, and docs/src/error_analysis.md for the full discussion + a worked example.

contour_df_samples(m::Minuit; kwargs...) -> DataFrame

DataFrame of the accepted Monte-Carlo parameter sets from get_contours_samples: one row per kept set, one column per free parameter (named after the parameters), plus a :delta_chisq column with each set's true Δχ². Requires using DataFrames (provided by the NativeMinuitDataFramesExt package extension).

Accepts all get_contours_samples keyword arguments.

contour_parameter_sets(ce::ContoursError) -> Vector{Vector{Float64}}

The full parameter set at every contour boundary point — one n-dimensional vector per point in ce.points, each holding the 2 contour coordinates (par_x, par_y) together with the n−2 other parameters at their profiled (re-minimized) values along the contour.

This is the native-Julia analogue of IMinuit.jl's get_contours, but with no extra fits: contour_exact already runs an inner re-minimization at each boundary point, so the full state is captured for free (IMinuit.jl re-fit each point only because PyCall could not return the inner MinimumState across the Python boundary).

Each vector is in the same coordinate frame as ce.points: physical (external) coordinates when the ContoursError came from the high-level contour(m, …) (which externalizes its output), and internal (sin/√-transformed) coordinates when it came from the low-level contour_exact(fmin, cf, …) on a bounded fit (identical to external for unbounded fits). Empty when ce came from the ellipse contour or when the contour was invalid.

Example

ce = contour_exact(fmin, cf, 1, 2; npoints = 24)
psets = contour_parameter_sets(ce)        # 24 full parameter vectors
# Re-evaluating the FCN at any set returns ≈ fmin + up (the boundary):
cf(psets[1]) ≈ fval(fmin) + cf.up

extremize(m::Minuit, f; cl=1, seeds=nothing, kwargs...) -> ExtremizeResult

Profile (Δχ²-region) confidence interval of a derived scalar f(θ):

[min f(θ), max f(θ)]   over   { θ : FCN(θ) ≤ m.fval + delta_chisq(cl, 1)·m.up }

with all free parameters varied simultaneously (fixed parameters stay pinned, limits are honoured). f receives the full EXTERNAL parameter vector, in m's parameter order. This is "MINOS for a function": for f(θ) = θ[i] it reproduces the MINOS interval, and in the linear-Gaussian limit it equals the projection-theorem result f̂ ± √(delta·cᵀCc) with the full parameter correlations included. Returns an ExtremizeResult; m is not mutated.

Why ndof = 1 even though all parameters move

The threshold is delta_chisq(cl, 1) because the quoted statement is ONE number: re-parametrize so f is itself a coordinate and this is a single-parameter interval with the others profiled out (Wilks: 1 dof). Using the joint delta_chisq(cl, n_free) here would over-cover (report errors up to ~3× too wide for 9 parameters). For a genuinely joint statement, pass an explicit delta — e.g. tracing a 2-D support function with delta = delta_chisq(cl, 2).

Multiple seeds are load-bearing, not an optimization

Each direction is an exterior-penalty minimization (obj = ∓f + λ·max(0, (FCN−bound)/up)², solved as a short penalty- continuation ladder of MIGRADs, λ = 1 → 100 → lambda warm-started) run from every seed; a single-seed run can stop short of the extremum and report a too-narrow interval in two distinct ways: (i) on a strongly correlated / ill-conditioned region the best-fit-anchored penalty stalls on the flat axis at a feasible but NON-extremal boundary point, and (ii) when the region splits into several disconnected low-χ² corridors the penalty cannot cross the barrier between them.

To cure (i) at no user effort, the result is, by default, floored/ceiled by the directional (HESSE-ellipse) endpoints θ̂ ± √δ·C∇f/σ_f (the :directional construction below): they are feasible and exact in the linear-Gaussian limit, so the reported interval is never narrower than the directional one, no matter what the penalty did (including a degenerate lambda that accepts nothing). It costs ONE extra directional probe (≈ n_free gradient + a dozen FCN + 2 f evaluations), NOT extra penalty seeds. Disable with directional_floor = false. If the direction is un-computable (∇fᵀC∇f ≤ 0, degenerate C, non-finite f at the probe) the floor is silently skipped. The floor does NOT cure case (ii): disconnected corridors are unreachable from a straight ray, so pass everything you have that touches other corridors (MCMC/ensemble members extreme in f, other find_solution_modes representatives) via seeds. Audit r.diagnostics: per-seed acceptance and f values, which seed won each side (winner_* == 0 means no penalty fit beat the best-fit value), and directional_floor.lo/.hi (whether the floor supplied that endpoint).

Cheap mode for expensive FCNs, and the f-failure contract

For an FCN/f that costs seconds, the default :full algorithm (multi-seed × ladder × rounds MIGRADs) can be hours per call — use mode = :directional (below) for the common near-linear case (≈ n_free + ~15 paired calls, ~50× cheaper), and on the :full path set rounds = 1, a small maxfcn, and strategy = 0. f may throw OR return a non-finite value at infeasible θ — both are safe: such probes become a finite high plateau the optimizer steers around (never NaN into MIGRAD), tallied as f_nonfinite in the diagnostics. A genuinely non-finite-f region may legitimately NARROW the interval (the optimizer cannot reach past it) — that is safe, not a bias. Do NOT instead return a sentinel like 0.0 from a failing f: that centers the endpoint (a silent bias). In mode = :directional, a non-finite f at a boundary crossing collapses that side onto the best fit, with a warning and a f_failed_lo/f_failed_hi diagnostic flag.

Keyword arguments

• cl::Real = 1 — confidence level, delta_chisq convention: cl ≥ 1 is nσ (1 → 68.27 %, 2 → 95.45 %), 0 < cl < 1 a probability (0.95 → 95 %). Threshold: Δχ² = delta_chisq(cl, 1).
• mode::Symbol = :full — :full is the multi-seed penalty extremization (handles non-linear / multi-corridor regions). :directional is the fast linear-direction crossing: it forms d = C·∇f at θ̂ (the Lagrange/projection direction this docstring's formula uses), secant/bisects the true FCN to its bound crossing on each side, and reports the true f there. Cost ≈ n_free (gradient) + ~2× a dozen FCN calls + 2 f calls — exact in the linear-Gaussian limit, and r.mode === :directional flags it so it is not mistaken for the full profile. It ignores seeds, does not chase non-linear corridors, and does not honour parameter limits that bind before the crossing — r.plo/r.phi may then lie OUTSIDE the limits (it warns when free parameters are bounded). The recommended workflow is :directional first, then :full (optionally seeded) only if you suspect non-linearity, a binding limit, or the two disagree. Requires m.covariance (run migrad!/hesse! first).
• grad_f = nothing — optional θ -> ∇f(θ) (full external length) for mode = :directional; replaces the n_free-call forward-difference gradient. Ignored by mode = :full.
• seeds = nothing — extra start points: a vector of full external parameter vectors, the rows of a matrix, or a single vector. The best fit is always prepended as seed 1. Fixed coordinates are re-pinned and free ones clamped into limits before fitting.
• directional_floor::Bool = true (mode = :full only) — floor/ceil the result by the directional (HESSE-ellipse) endpoints so it is never narrower than the linear-Gaussian interval (see "Multiple seeds are load-bearing"). Uses a numerical gradient (ignores grad_f, like the rest of :full); set false to skip the extra directional probe.
• lambda::Real = 1e4 — FINAL penalty stiffness; the continuation ladder is unique([min(lambda,1), min(lambda,100), lambda]). The raw optimum overshoots the boundary by O(1/lambda) and is then pulled back onto it.
• accept_tol::Real = 0.05 — acceptance gate, in units of up: a penalty optimum with FCN > bound + accept_tol·up is discarded as not converged onto the region. The gate is applied to the raw penalty optimum, BEFORE the boundary pull-back — and that optimum always overshoots the boundary by O(1/lambda), so accept_tol must stay above the overshoot (rule of thumb: ≳ 10/lambda in up units). In particular accept_tol = 0 is essentially never satisfiable: it rejects every candidate and silently collapses the result to the best-fit value — tighten lambda (stiffer penalty ⇒ smaller overshoot), not the gate. Exact feasibility of the reported endpoints is the pull-back's job, not the gate's; certify it via the fcn_* diagnostics.
• delta::Union{Real,Nothing} = nothing — explicit Δχ² threshold override (FCN units of up); when given, cl is ignored (and recorded as NaN).
• rounds::Integer = 4 — maximum warm-started repeats of the penalty ladder per seed and direction. Inside the region the objective is just ∓f (zero penalty, near-zero curvature for a near-linear f), so one chain can stall part-way through a long interior traverse — e.g. along a curved χ² valley — at a feasible but non-extremal point; repeats advance it. Iteration stops as soon as a round no longer improves f (a smooth problem therefore runs 2 rounds: one to converge, one to confirm).
• strategy = m.strategy, maxfcn = nothing — per-penalty-fit MIGRAD strategy and call budget. For an expensive FCN, strategy = 0 (no inner HESSE) and a modest maxfcn cut the cost sharply.
• iterate::Integer = 5 — the per-penalty-MIGRAD retry budget forwarded to migrad! (its _robust_low_level_fit parity). The default helps a stiff-penalty stage that stalls invalid recover; for an expensive FCN set iterate = 1 to forbid retries (the cheapest setting — combine with rounds = 1).
• on_unit = nothing — a callback u -> … fired once per completed penalty- MIGRAD unit with u = (side, seed, round, stage, cur, nfcn, fcn, valid) — the granularity of one MIGRAD (≈ one expensive step). Use it for live progress or to checkpoint partial work externally (a kill then only loses the in-flight unit). Attaching it adds one FCN evaluation per unit (to report fcn); the default nothing pays nothing. Ignored by mode = :directional.

Cost & guarantees

2 directions × (number of distinct seeds) × (ladder stages, ≤ 3) × (rounds actually run — ≥ 2 whenever rounds ≥ 2: one to converge plus one to confirm; exactly 1 under rounds = 1, the profile_band default) MIGRAD runs of the penalty objective (each evaluation calls the FCN and f once), plus ≤ ~50 FCN calls per accepted endpoint for the boundary pull-back. lo ≤ f(θ̂) ≤ hi by construction; endpoints satisfy FCN ≤ bound exactly when the pull-back applies (the typical case) and FCN ≤ bound + accept_tol·up always. A side on which no penalty fit is accepted falls back to the best-fit value with a warning — treat that side as failed and investigate the diagnostics.

Example

m = Minuit(chi2, x0; names = names, limits = limits)
migrad!(m)

# 68.3 % interval for a derived quantity (here: the model curve at x = 15)
r = extremize(m, θ -> θ[1] + θ[2] * 15.0)
r.lo, r.hi          # the interval
r.plo, r.phi        # parameter vectors realizing the endpoints
r.diagnostics       # per-seed audit: who converged, who was accepted, who won

# 95 % (2σ), seeding from ensemble members extreme in f
r2 = extremize(m, f; cl = 2, seeds = ens[sortperm(f.(eachrow(ens)))[[1, end]], :])

See also profile_band (pointwise band of a curve family), minos! (the f(θ) = θᵢ special case), delta_chisq, and docs/src/error_analysis.md for the full decision guide.

profile_band(m::Minuit, f, xs; cl=1, seeds=nothing, warm=true, passes=2,
             include_best=true, kwargs...) -> ProfileBand

Pointwise profile-likelihood error band of a curve family f(x, θ) on the grid xs — x first, θ the full EXTERNAL parameter vector in m's parameter order: the package-wide model(x, …) convention, and the same callback shape as quantile_band. At each grid point, extremizes θ -> f(x, θ) over the same fixed region {FCN ≤ m.fval + delta_chisq(cl, 1)·m.up} (all free parameters varied, limits/fixed honoured). Returns a ProfileBand with the envelope (lo, hi), the per-point extremal parameter vectors (plo, phi), the best-fit curve (fbest), the failure count and per-point diagnostics. m is not mutated.

This is the standard pointwise (Δχ² ≤ 1 at cl = 1) construction for figure bands: each x carries its own cl statement, and the band contains the best-fit curve by construction — unlike posterior-quantile bands, which can exclude it when a parameter sits on a limit.

Sweep strategy

The grid is swept passes times, alternating forward/reverse, keeping the better (more extreme, still feasible) envelope per point. With warm = true each side carries the previous point's extremal parameter vector as an extra seed — extremal points move continuously along a smooth curve family, so the warm seed is usually the best one. The full seed pool (best fit + seeds) is also used at every point: corridor coverage must not depend on the sweep having visited the right corridor earlier. For an expensive FCN, the cost model is

#MIGRAD ≈ length(xs) × 2 sides × passes × (pool + warm + incumbent) × stages × rounds

(stages ≤ 3, the penalty-continuation ladder; rounds defaults to 1 here) — trim seeds, lower passes, or set maxfcn (a per-MIGRAD budget) to control it.

Keyword arguments

cl, seeds, lambda, accept_tol, delta, strategy, maxfcn as in extremize, plus:

• mode::Symbol = :full — :full runs the multi-seed penalty extremization at every grid point; :directional runs the fast C·∇f crossing at every point (see extremize's mode). For an expensive FCN and a near-linear curve family, :directional makes the whole sweep orders-of-magnitude cheaper (≈ npoints × (n_free + ~15) evaluations) and is exact in the linear-Gaussian limit; it ignores seeds/warm/passes/ rounds/iterate/on_unit, does not honour limits, and flags a point that fails (non-finite f or un-computable direction) in nfail/diagnostics with a best-fit fallback. b.mode records which was used. Requires m.covariance.
• grad_f = nothing — optional (x, θ) -> ∇_θ f(x, θ) (full external length) for mode = :directional; replaces the per-point forward-difference gradient. Ignored by mode = :full.
• rounds::Integer = 1 — per-point penalty-ladder repeats (see extremize). The band default is 1 because the sweep itself iterates each point: the warm seed, the stored incumbent and the forward/reverse passes re-polish every point several times. Raise it if isolated points lag behind their neighbours.
• warm::Bool = true — warm-start each point from the neighbouring point's extremal parameters (per side, per pass direction).
• passes::Integer = 2 — sweep count (1 = forward only; 2 adds the reverse sweep; more keep alternating). Two passes let a corridor discovered mid-grid propagate to BOTH sides.
• include_best::Bool = true — keep the best-fit value as a zero-cost feasible candidate per point, guaranteeing lo ≤ fbest ≤ hi and a finite band even where every penalty fit failed (such failures are still counted in nfail/diagnostics). With false, a fully-failed side is NaN — useful when benchmarking the optimizer itself.
• iterate::Integer = 5 — per-penalty-MIGRAD retry budget (see extremize); iterate = 1 is the cheapest for an expensive FCN.
• on_unit = nothing — per-MIGRAD-unit progress callback (see extremize); the record additionally carries (x, point, pass) for the grid point. Enables external checkpointing of a long band sweep.
• verbose::Bool = false — @info one line per pass.

Example

m = Minuit(chi2, x0; names = names); migrad!(m)
mgrid = 4360.0:2.0:4520.0
band  = profile_band(m, (x, θ) -> moment_P2(x, θ), mgrid;
                     seeds = ens_extremes)         # ensemble extreme members
band.nfail == 0 || @warn "inspect band.diagnostics"
# plot: fill between band.lo and band.hi, line at band.fbest

See also extremize, mnprofile (profile of a single parameter), and docs/src/error_analysis.md.

ExtremizeResult

Result of extremize: the profile (Δχ²-region) interval of a derived scalar f(θ), the extremal parameter vectors realizing it, and the per-seed audit trail.

Fields

• lo::Float64, hi::Float64 — the interval endpoints [min f, max f] over the region FCN ≤ bound. In the typical case each is attained at a FEASIBLE point (FCN ≤ bound exactly, after the boundary pull-back); a candidate whose pull-back is not locally possible is kept raw, within FCN ≤ bound + accept_tol·up (its record has projected = false).
• plo, phi — full external parameter vectors at which lo/hi are attained (fixed parameters pinned at their values).
• fbest::Float64 — f at the best fit; lo ≤ fbest ≤ hi by construction.
• bound::Float64 — the FCN acceptance bound m.fval + delta·up.
• delta::Float64 — the Δχ² threshold actually used (delta_chisq(cl, 1), or the explicit delta override).
• cl::Float64 — the confidence-level argument as given (NaN when an explicit delta override was used).
• up::Float64 — the fit's error definition (m.up).
• mode::Symbol — :full (the multi-seed penalty extremization) or :directional (the fast linear-direction crossing; see extremize). The diagnostics schema differs between the two modes.
• diagnostics::NamedTuple — for mode = :full: (min, max, winner_min, winner_max, naccepted_min, naccepted_max, fcn_min, fcn_max, directional_floor). min/max are per-seed record vectors, one row per penalty fit, with fields seed (index into the seed pool; seed 1 is always the best fit), converged (MIGRAD validity), accepted (passed the bound + accept_tol·up gate), projected (was pulled back onto the boundary), fcn (FCN at the raw penalty optimum), f_raw/f (f before/after the pull-back), nfcn, and f_nonfinite (count of probes where f threw or returned non-finite and was steered around — see the f-failure contract under extremize). winner_min/winner_max give the seed index whose candidate won each side — 0 means no penalty fit beat the best-fit value itself, which has two very different readings, disambiguated by naccepted_*: with naccepted_* > 0 the best fit is genuinely extremal for that side (e.g. f's unconstrained optimum sits at θ̂ — healthy); with naccepted_* == 0 every penalty fit was rejected and that side FAILED (a warning fired). fcn_min/fcn_max are the FCN values at plo/phi: compare with bound to certify each endpoint's feasibility directly. Inspect these to audit seed coverage (the multi-corridor failure mode described in the docstring). directional_floor is (lo::Bool, hi::Bool) — whether the directional floor/ceiling (not the penalty) supplied that endpoint. For mode = :directional: (grad, dir, gCg, alpha_lin, alpha_lo, alpha_hi, fcn_lo, fcn_hi, nfcn, nf) — the f-gradient and search direction C·∇f at θ̂, the linear step √(delta/gᵀCg), the two true-FCN boundary crossings alpha_lo/alpha_hi (± along dir) with their FCN values (≤ bound), and the FCN/f call counts.

ProfileBand

Result of profile_band: the pointwise profile-likelihood envelope of a curve family f(x, θ) on a grid.

Fields

• x::Vector{Float64} — the grid (a copy of xs).
• lo, hi — the band edges: per point, [min f, max f] over the region FCN ≤ bound. As in extremize, an edge typically sits at a feasible point (FCN ≤ bound exactly, after the boundary pull-back) and always within FCN ≤ bound + accept_tol·up — check the per-point fcn_lo/fcn_hi diagnostics to certify. NaN only when include_best = false and every penalty fit at that point/side was rejected.
• plo, phi — per-point extremal parameter vectors (nothing exactly when the corresponding edge is NaN).
• fbest::Vector{Float64} — the best-fit curve f(x, θ̂); with include_best = true (default) lo .≤ fbest .≤ hi by construction.
• bound, delta, cl, up — as in ExtremizeResult.
• mode::Symbol — :full (per-point multi-seed penalty extremization) or :directional (per-point fast C·∇f crossing; see profile_band). The per-point diagnostics schema differs between the two modes.
• nfail::Int — number of failed (point, side) groups: for :full, (point, side, pass) extremization groups in which NO penalty fit passed the acceptance gate; for :directional, (point, side) crossings where f was non-finite or the direction was un-computable. The best-fit fallback still keeps the band finite when include_best = true; 0 on a healthy sweep.
• diagnostics — per-point NamedTuples. :full: (x, failed_lo, failed_hi, accepted_lo, accepted_hi, nfits_lo, nfits_hi, fcn_lo, fcn_hi) — cumulative accepted / attempted penalty-fit counts over all passes, failed_* flags for a side that NEVER had an accepted fit, and the FCN values at the stored extremal points. :directional: (x, failed_lo, failed_hi, fcn_lo, fcn_hi, gCg, nfcn, nf) — the per-point feasibility flags, boundary FCNs, ∇fᵀC∇f, and FCN/f call counts.

mcmc_sample(m::Minuit; nsteps=52_000, burn=2_000, thin=25,
            proposal=:hesse, scale=0.3, target_accept=nothing,
            adapt_every=100, seed=nothing, rng=nothing, warn=true)
    -> LikelihoodEnsemble

Sample parameter sets from the fit's likelihood L(θ) ∝ exp(−fcn(θ) / (2·up)) with a random-walk Metropolis chain on the exact FCN — χ²: exp(−Δχ²/2); negative log-likelihood (up = 0.5): exp(−Δ(−log L)) — started at the current best fit. The returned LikelihoodEnsemble feeds quantiles / quantile_band (marginal quantile intervals and pointwise bands for any derived quantity) and can be stored with save_ensemble as a reusable error set.

This is the likelihood-ensemble leg of the error-triangulation (profile/MINOS extremization ↔ ensemble quantiles); iminuit has no native analogue (Python users attach emcee). Call migrad!(m) (and ideally hesse!(m), for a well-shaped proposal) first.

NOT the same as get_contours_samples

get_contours_samples samples the confidence region Δχ² ≤ delta_chisq(cl, ndof) — a hard cut, giving region extents. This function runs a posterior chain: no Δχ² cut at all; samples concentrate where the likelihood mass is, i.e. at Δχ² ≈ n_free (the high-dimensional volume effect — for 9 free parameters P(Δχ² ≤ 1) ≈ 5.6e-4, so a likelihood chain essentially never visits the Δχ² ≤ 1 region, and does not need to). Use the region sampler for joint confidence regions; use this chain for likelihood-weighted quantiles of derived quantities.

Marginal quantiles vs profile envelopes — read before quoting bands

A (16%, 84%) quantile band built from this ensemble is a marginal (posterior-mass) construction. The profile construction (mnprofile, MINOS, or constrained extremization of a derived quantity over Δχ² ≤ 1) contains the best fit by construction; a marginal quantile band need not. The two agree in the near-Gaussian interior and separate legitimately at parameter limits: the truncation piles the posterior mass one-sidedly, so the band can shift away from a best fit sitting on the boundary (mode ≠ median). That separation is a property of the constructions, not a sampler failure, and more samples will not make it go away. See docs/src/error_analysis.md for the full comparison table.

Algorithm

Plain Metropolis: propose θ' = θ + scale·δ with symmetric δ, accept with probability min(1, exp(−(fcn(θ′) − fcn(θ)) / (2·up))), record every thin-th state after the first burn steps (nkept = (nsteps − burn) ÷ thin). Fixed parameters never move. Proposals outside parameter limits are rejected before the FCN is called — the chain samples the likelihood truncated to the allowed box (see above). Non-finite FCN values are likewise rejected, never accepted. The chain bypasses the fit's call counter (m.nfcn is untouched) and never mutates m.

The proposal shape only affects mixing efficiency, never what the chain converges to (any symmetric proposal has the same stationary distribution) — so an imperfect Σ is far less dangerous here than in get_contours_samples, where proposal under-coverage biases the region.

Keyword arguments

• nsteps=52_000, burn=2_000, thin=25 — chain length, discarded burn-in, and keep-every-thin stride; the defaults yield 2000 kept sets (a field-validated recipe). Thinning mainly buys decorrelation; for quantile work ~1–2 k kept sets are usually plenty.
• proposal=:hesse — the proposal step shape, one of
  • :hesse — multivariate Gaussian with the fit covariance (correlations included). Falls back to :errors with a warning when no covariance is available or it looks unreliable (invalid fit / forced pos-def / failed HESSE);
  • :errors — independent per-coordinate Gaussians with σ = the fit's parabolic errors as frozen in the fit result (m.fmin.ext_errors — numerically m.errors right after the fit) — the classic hand-rolled choice;
  • an AbstractVector of length n_free — explicit per-free-parameter σ (external units). The escape hatch when both the covariance and the parabolic errors are meaningless, e.g. a parameter pinned at a limit;
  • an AbstractMatrix (n_free × n_free) — explicit proposal covariance.
• scale=0.3 — overall step multiplier: step σ = scale × (proposal σ). Field experience: 0.25–0.35 × HESSE σ gives acceptance ≈ 0.2–0.3 for ~10 free parameters. (Random-walk theory suggests ≈ 2.38/√n_free for a perfectly matched :hesse proposal; in practice start at 0.3 and let target_accept tune. In LOW dimension the same scale accepts much more — ≈ 0.8 for a 3-parameter Gaussian fit — which is conservative, not a misconfiguration; thinning absorbs the extra correlation.)
• target_accept=nothing — when set (e.g. 0.25), the scale is adapted toward that acceptance during burn-in only (every adapt_every steps, multiplicative update clamped to ×[0.5, 2] per round), then frozen, so the kept chain has a fixed kernel. Requires burn ≥ adapt_every.
• seed — make the chain reproducible (seed=11 etc.); rng — supply your own AbstractRNG instead (mutually exclusive with seed).
• warn=true — emit tuning warnings (unreliable-covariance fallback, post-burn acceptance < 0.05 or > 0.9).

Returns

A LikelihoodEnsemble: kept sets in .samples (full external vectors, fixed columns constant), FCN values in .fvals, post-burn .acceptance, the final .scale, and the chain metadata.

Example

m = Minuit(chi2, x0; names = names, limits = limits)
migrad!(m); hesse!(m)

ens = mcmc_sample(m; seed = 11)                    # 52 k steps → 2000 sets
ens.acceptance                                     # aim for ≈ 0.2–0.3

# scalar derived quantity: 16/50/84% quantiles
q16, q50, q84 = quantiles(ens, θ -> θ[1] - θ[2])

# pointwise 16–84% band of a curve
band = quantile_band(ens, (x, θ) -> model(x, θ), xgrid)

save_ensemble("ensemble.dat", ens; comment = "error set B")  # reusable

See also quantiles, quantile_band, save_ensemble, load_ensemble, get_contours_samples, and docs/src/error_analysis.md.

LikelihoodEnsemble

Likelihood-weighted parameter ensemble produced by mcmc_sample (or reloaded by load_ensemble). One row per kept Metropolis step, in full external coordinates (fixed parameters appear as constant columns), so any user model f(θ_full) evaluates directly on a row.

Fields

• samples::Matrix{Float64} — nkept × npar kept parameter sets (row = one set, columns ordered as m.parameters).
• fvals::Vector{Float64} — the FCN value (χ² or negative log-likelihood, whatever the fit minimized) at each kept set. Per-sample χ²-equivalent displacements are (fvals .- fbest) ./ up.
• names::Vector{String} — all parameter names (length npar).
• free::Vector{Bool} — which columns were varied by the chain.
• best::Vector{Float64} — the chain's start point (the fit's best values, full external vector).
• fbest::Float64 — the FCN evaluated at best when the chain started.
• up::Float64 — the fit's errordef (1 for χ², 0.5 for −log L).
• acceptance::Float64 — accepted fraction of the post-burn-in proposals (healthy random-walk Metropolis: ≈ 0.2–0.4).
• nsteps::Int, burn::Int, thin::Int — the chain settings; the kept rows are (nsteps − burn) ÷ thin per chain, so the total is that times the number of chains (or walkers, for the :stretch ensemble).
• scale::Float64 — the final proposal scale (after any burn-in adaptation; equals the input scale when target_accept was not used). NaN for the :nuts sampler, which has no proposal scale.
• proposal::Symbol — the proposal actually used: :hesse, :errors, :steps (explicit per-coordinate σ) or :matrix (explicit covariance); :unknown for ensembles loaded from a foreign file.
• seed::Union{Nothing,UInt64} — the RNG seed, when one was given.

Behaves like a collection of parameter vectors

length(ens) is the number of kept sets; ens[i] returns the i-th parameter vector (a copy); iteration yields the rows, so [f(θ) for θ in ens] evaluates a derived quantity over the whole ensemble (which is exactly what quantiles and quantile_band automate).

A chain can wander below the fit minimum (multi-modal landscape, or an unconverged fit): minimum(ens.fvals) < ens.fbest is then a finding, not an error — re-minimize (see find_deeper_minimum) before quoting errors about the old minimum. The show method flags this.

quantiles(ens::LikelihoodEnsemble, f; p=(0.16, 0.5, 0.84), warn=true)
    -> Vector{Float64}

Quantiles of a scalar derived quantity f(θ) over the likelihood ensemble: evaluates f on every kept parameter set (full external vector) and returns the quantiles at probabilities p, in order.

The default p gives the 16% / median / 84% triplet — a marginal "1σ-equivalent" interval. Include 0 / 1 in p for the ensemble minimum / maximum. Non-finite f values are dropped (with a warning).

These are marginal (posterior-mass) quantiles — see the mcmc_sample docstring for how they relate to (and legitimately differ from) profile/MINOS intervals, especially at parameter limits.

q16, q50, q84 = quantiles(ens, θ -> θ[2] - θ[1])
lo, hi        = quantiles(ens, θ -> θ[1]; p = (0.16, 0.84))

See also quantile_band for pointwise bands of a curve.

quantile_band(ens::LikelihoodEnsemble, f, xs; p=(0.16, 0.84),
              curve=false, warn=true) -> Matrix{Float64}

Pointwise quantile band of a curve over the likelihood ensemble: at each grid point xs[i], the quantiles (at probabilities p) of f(xs[i], θ) across all kept parameter sets θ. Returns a length(xs) × length(p) matrix — with the default p, column 1 is the 16% (lower) and column 2 the 84% (upper) band edge.

• f — f(x, θ)::Real with θ the full external parameter vector. With curve = true, instead f(θ)::AbstractVector returning the whole curve over xs in one call (one f call per ensemble member — use this when evaluating the model pointwise repeats expensive shared work).
• p — quantile probabilities; e.g. p = (0.025, 0.16, 0.5, 0.84, 0.975) for median + 1σ + 2σ-equivalent bands.
• Non-finite values are dropped per grid point (with a warning).

This is a marginal quantile band (likelihood-mass construction), the companion of — not a substitute for — the profile envelope band (pointwise extremization of f over Δχ² ≤ delta_chisq(cl, ndof)): the profile band contains the best-fit curve by construction, a quantile band need not (boundary effects shift the mass one-sidedly; mode ≠ median). Quote which construction you used. See mcmc_sample and docs/src/error_analysis.md.

band = quantile_band(ens, (x, θ) -> model(x, θ), xgrid)
plot(xgrid, mid; ribbon = (mid .- band[:, 1], band[:, 2] .- mid))

# expensive model: one call per ensemble member returns the whole curve
band = quantile_band(ens, θ -> model_curve(xgrid, θ), xgrid; curve = true)

save_ensemble(path_or_io, ens::LikelihoodEnsemble; comment="")

Write the ensemble as plain text: #-comment header (metadata + optional user comment), then one line per kept set — the FCN value followed by all parameter values, space-separated. Floats are written in shortest round-trip form, so load_ensemble reproduces them exactly.

The format is deliberately the time-honoured hand-rolled one (# header

• fval p₁ p₂ … rows): existing files of that shape load fine, and the

saved file is directly consumable by gnuplot / numpy / a 5-line Julia loop. Ensembles are worth saving: any future derived quantity gets its band by evaluating over the stored sets — no re-sampling.

Because the columns are whitespace-separated, parameter names must not contain whitespace (they are written to the # names: header); a name with a space throws rather than corrupting the round-trip.

save_ensemble("ensemble_B.dat", ens; comment = "error set B, ρ = -0.295")

load_ensemble(path_or_io; names=nothing, up=nothing) -> LikelihoodEnsemble

Read an ensemble written by save_ensemble — or any plain-text file of the same shape: # comment lines, then one row per sample, fval p₁ p₂ … space-separated (the classic hand-rolled ensemble format).

Metadata found in # key: value header lines (as written by save_ensemble) is restored; anything missing gets a placeholder (names = ["p1", …], free all true, best/fbest/up/… NaN, proposal = :unknown). Override with the names / up keywords when loading a foreign file. Quantile analysis (quantiles / quantile_band) only needs the samples, so it works on foreign files as-is.

bayesian(m::Minuit; prior=:flat, level=0.6827, interval=:central, kwargs...)
    -> BayesianReport

One-step, non-mutating Bayesian analysis convenience wrapper. It calls posterior_sample, summarizes the result, and returns a BayesianReport. It never writes Bayesian intervals into m.errors, m.covariance, or MINOS state.

BayesianReport

Non-mutating one-step Bayesian report returned by bayesian. It holds the reusable PosteriorSample, the requested credibility level, the interval method, and the parameter summary table.

posterior_sample(m::Minuit; prior=:flat, kwargs...) -> PosteriorSample
posterior_sample(prob::PosteriorProblem; kwargs...) -> PosteriorSample

Sample the posterior exp(-fcn/(2*up)) * prior in external coordinates. The returned object is independent of m; the fit state and m.nfcn are not mutated. Three samplers are available. The default, :stretch, is gradient-free and affine-invariant — it works for any FCN (including a non-differentiable complex-buffer χ²) and mixes well on strongly correlated posteriors:

• sampler = :metropolis — a random-walk Metropolis chain, nchains of them (default 4). The proposal is set by proposal (:hesse / :errors / a σ vector / a covariance), tuned by scale and optional target_accept. Chains 2…n start over-dispersed at overdisperse × the proposal/HESSE scale from the best fit (default 2, ≈2σ wider than the posterior) — the condition that makes the multi-chain split-R̂ a real convergence test.
• sampler = :stretch (the default) — the affine-invariant ensemble sampler (Goodman & Weare; the emcee kernel). nwalkers walkers (default max(2·n_free+2, 8), rounded up to even) explore the posterior with stretch moves of scale stretch (a, default 2). It is gradient-free (works for any FCN, including ones that cannot be auto-differentiated) and handles strongly correlated / skewed posteriors far better than a single random-walk chain. For :stretch, nsteps counts ensemble iterations (each updates every walker; defaults nsteps=6000, burn=1000, thin=10), and each walker is treated as a chain for the R̂ / ESS diagnostics. proposal / scale / target_accept / overdisperse are random-walk-only and ignored here.
• sampler = :nuts — gradient-based NUTS (No-U-Turn HMC), provided by the AdvancedHMC extension (load AdvancedHMC, LogDensityProblems, LogDensityProblemsAD, TransformVariables, ForwardDiff alongside NativeMinuit). Bounded parameters are mapped to unconstrained ℝ with the proper log-Jacobian and sampled with a ForwardDiff gradient; it is the most efficient sampler for smooth, higher-dimensional posteriors but requires an auto-differentiable FCN (no finite-difference fallback — it errors and points to :stretch) and cannot start from a best fit on a parameter limit. All built-in cost objects (LeastSquares / UnbinnedNLL / ExtendedUnbinnedNLL / BinnedNLL / ExtendedBinnedNLL / CostSum) are ForwardDiff-differentiable and work with :nuts, provided the model / pdf / cdf you supply is itself auto-differentiable — a user function that hard-codes Float64 internally (e.g. a complex-buffer χ²) still needs the gradient-free :stretch. Defaults nsteps=2000 (per chain, incl. warmup), burn=1000, target_accept=0.8.

The reported rhat is the basic split-R̂ (not rank-normalized / folded); for skewed or boundary-truncated marginals also check effective_sample_size and the trace rather than trusting R̂ < 1.01 alone. Bayesian credibility levels are plain probabilities; no frequentist nσ overload is used.

PosteriorProblem(m::Minuit; prior=:flat)

Snapshot of a Minuit fit plus an explicit prior, defining the posterior target

log p(θ | data) = -fcn(θ) / (2 * up) + logprior(θ)

in full external parameter coordinates. The object is non-mutating and stores a snapshot of the fit state; use isconsistent to compare it with a later Minuit object.

MVP semantics: Minuit limits are physical posterior support, intersected with the prior support. If the current best point is outside that effective support, construction fails loudly rather than starting a dead chain.

PosteriorSample

Posterior sample returned by posterior_sample. The ensemble field is a regular LikelihoodEnsemble in full external coordinates; its fvals remain likelihood FCN values. Posterior-specific provenance lives in prior, loglik_kept, logpost_kept, chain IDs, diagnostics, and warnings.

isconsistent(prob::PosteriorProblem, m::Minuit) -> Bool

Return whether m still matches the fit snapshot captured in prob.

Prior

Log-prior provenance for posterior_sample / bayesian. logdensity(θ_full) receives the full external parameter vector and returns an unnormalized log density. support_lo/support_hi are full-length support vectors (±Inf for unbounded coordinates). informative marks parameters with explicit non-flat prior components and is used by combine_priors to reject accidental overlap.

all_components_proper is provenance only; it is not a proof that the posterior is proper.

flat_prior(m::Minuit) -> Prior

Flat prior in NativeMinuit's full external parameter coordinates. This is not "no prior"; it is a parameterization-dependent coordinate choice. On unbounded coordinates posterior propriety relies on the likelihood.

normal_prior(m::Minuit, par, μ, σ) -> Prior

Gaussian prior on parameter par, flat on the remaining coordinates. The log density is unnormalized: -0.5*((x-μ)/σ)^2.

uniform_prior(m::Minuit, par, lo, hi) -> Prior

Proper uniform prior support for one parameter, flat on the remaining coordinates. Returns -Inf outside [lo, hi].

half_normal_prior(m::Minuit, par, σ) -> Prior

Half-normal prior for a non-negative or lower-bounded parameter. If par has a finite Minuit lower limit, the half-normal is centered at that lower limit; otherwise it is centered at zero and supports par >= 0.

combine_priors(p1, p2, ...) -> Prior

Combine disjoint informative prior components by adding their log densities and intersecting their supports. MVP behavior is deliberately strict: two informative components on the same parameter raise an error.

credible_interval(post, par; level=0.6827, method=:central)

Equal-tailed marginal Bayesian credible interval for parameter par. MVP supports method=:central only; HPD intervals are deferred because multimodal marginals may have disjoint highest-density regions.

derived_interval(post, f; level=0.6827, method=:central, warn=true)

Sample-wise credible interval for a scalar derived quantity f(θ_full). No delta method or linear propagation is used. warn=false silences the "dropped N non-finite values" notice from the underlying quantile pass.

upper_limit(post, par; level=0.90) -> CredibleLimit

One-sided Bayesian upper credible limit: the value x with P(par <= x | data, prior) = level.

lower_limit(post, par; level=0.90) -> CredibleLimit

One-sided Bayesian lower credible limit: the value x with P(par >= x | data, prior) = level.

CredibleLimit

One-sided Bayesian credible limit returned by upper_limit or lower_limit. This is not a CLs, Feldman-Cousins, MINOS, or profile confidence limit; it is conditional on the stated prior.

posterior_summary(post; level=0.6827)

Return per-parameter summary rows (parameter, mean, median, std, lower, upper, level) using central marginal credible intervals.

posterior_mean(post::PosteriorSample, par) -> Float64

Posterior mean of parameter par (name or index) over the kept samples.

posterior_median(post::PosteriorSample, par) -> Float64

Posterior median of parameter par (name or index).

posterior_std(post::PosteriorSample, par) -> Float64

Posterior standard deviation of parameter par (sample std, N − 1 divisor). This is a marginal posterior spread, not a HESSE/MINOS error. Returns NaN for ≤ 1 kept samples, and 0.0 for a constant column (e.g. a fixed parameter).

effective_sample_size(post::PosteriorSample, par) -> Float64

Effective sample size for parameter par, summed over chains (autocorrelation- adjusted). NaN for a fixed parameter. A small ESS relative to the kept-sample count means the chain mixed slowly; raise nsteps/thin or retune the proposal.

rhat(post::PosteriorSample, par) -> Float64

Split-R̂ convergence diagnostic for parameter par (needs nchains ≥ 2). Values near 1 (conventionally < 1.01) indicate the chains have mixed; NaN for a fixed parameter or a single chain. This is the basic split-R̂ (not rank-normalized / folded), so for a skewed or boundary-truncated marginal pair it with effective_sample_size and a trace check.

bootstrap(model::Function, data::Data, start; kws...) -> BootstrapResult

Bootstrap error analysis for a χ² fit of model(x, par) to data. start is either the initial-value AbstractVector (the full-data fit is run internally to obtain the anchor optimum θ̂_full) or a Minuit template (its configuration — names, limits, fixed flags, errordef, strategy, tol, precision — is cloned, and its current values seed the anchor; it is fitted first if not already). A user analytic grad is deliberately not carried into the re-fits (it is a closure over the original data, so it would be wrong for a resampled set); each resample uses the numerical gradient.

Keyword arguments:

• nresample::Int = 1000 — number of bootstrap resamples.
• kind::Symbol = :nonparametric — :nonparametric resamples the N data points with replacement and re-fits; :parametric regenerates yᵢ* = model(xᵢ, θ̂_full) + σᵢ·zᵢ, zᵢ ~ 𝒩(0,1), holding x and err fixed (Gaussian error model from data.err).
• seed::Union{Integer,Nothing} = nothing — master RNG seed; pass an integer for reproducible (and threaded == serial bit-identical) results. When nothing, rng is used as-is.
• rng::AbstractRNG = Random.default_rng() — master RNG when seed is not set.
• warm_start::Bool = true — seed each re-fit from θ̂_full (fast, recommended). false cold-starts every re-fit from the original start values.
• ci_level::Real = 0.68 — percentile-CI coverage (default ≈ ±1σ).
• covariance::Bool = false — also return the bootstrap covariance matrix.
• filter_invalid::Bool = true — compute summary stats over converged re-fits only. false includes every re-fit that produced a finite estimate (converged or not); a thrown / non-finite re-fit is dropped either way (it would poison the statistics). All rows, valid or not, are always kept in samples.
• threaded::Bool = false — run the re-fits across threads (requires julia -t N); the model closure must be thread-safe (pure). Deterministic: threaded and serial runs are bit-identical for the same seed.
• further kws... flow to model_fit / Minuit (e.g. name, limits) when start is a plain vector.

The returned BootstrapResult carries the θ̂ sample matrix plus per-parameter mean / std / percentile CIs. The bootstrap std ≈ HESSE error for a well-specified Gaussian fit and diverges when the error model is wrong.

bootstrap(refit::Function, data; kws...) -> BootstrapResult

Generic (interface-iii) nonparametric bootstrap. data need only support length(data) and data[idxvec] (e.g. Data, or any indexable collection), and refit(subdata) -> θ̂::AbstractVector re-fits a resampled dataset and returns the parameter vector. Warm-starting is the caller's responsibility (let the closure start from the full-data optimum).

Keywords: nresample, seed, rng, ci_level, covariance, threaded as for the model-based method, plus names::Union{Vector{String},Nothing}=nothing to label the parameters (default p1, p2, …). kind is fixed to :nonparametric here — a parametric bootstrap needs the model + error model, so use the model/Data method for that.

The initial refit(data) probe (used to size the parameter vector) is not guarded — if it throws, the whole call fails; only the per-resample re-fits are caught and recorded as invalid.

bootstrap(cost::AbstractCost, start; kws...) -> BootstrapResult

Bootstrap for a Julia-native cost object (interface i) — the cost carries its data + model/pdf, so no separate model is passed. Supported for the point-level costs LeastSquares, UnbinnedNLL, and ExtendedUnbinnedNLL (nonparametric resampling of the data points with replacement); LeastSquares additionally supports kind = :parametric (Gaussian y-regeneration from the best-fit model). Binned costs (BinnedNLL / ExtendedBinnedNLL) and the composite CostSum are not point-resamplable and raise an ArgumentError pointing to the generic interface.

start is an initial-value vector or a fitted Minuit; all keywords (nresample, kind, seed, rng, warm_start, ci_level, covariance, filter_invalid, threaded) and the returned BootstrapResult match the model + Data method. Masked points are excluded (the cost is resampled over its active set).

jackknife(model::Function, data::Data, start; kws...) -> JackknifeResult

Delete-1 (or delete-d block) jackknife for a χ² fit of model(x, par) to data. start is an initial-value vector or a Minuit template (configuration cloned), exactly as for bootstrap.

Keyword arguments:

• d::Int = 1 — block size. d = 1 deletes one point per re-fit (g = N re-fits); d > 1 deletes consecutive blocks (g = ceil(N/d) re-fits), preserving any serial correlation within a deleted block. The block (delete-d) jackknife is intended for serially-correlated data and is a deliberately COARSE estimator: with only g groups it has ≈ √(2/(g-1)) relative scatter, and for IID data consecutive blocks of sorted data are non-exchangeable (a mild downward variance bias) — shuffle the data before blocking if it is IID. The (g-1)/g variance scaling is exact for equal blocks (when d ∣ N) and approximate for the shorter trailing block otherwise. Prefer d = 1 unless the data are correlated.
• warm_start::Bool = true — seed each re-fit from θ̂_full.
• threaded::Bool = false — run the re-fits across threads (deterministic; the jackknife uses no randomness, so threaded and serial results are identical).
• further kws... flow to model_fit / Minuit when start is a vector.

The returned JackknifeResult reports the bias estimate (g-1)·(θ̄ - θ̂_full), the bias-corrected estimate, the jackknife variance ((g-1)/g)·Σ(θ̂₍ⱼ₎ - θ̄)², and the full jackknife covariance matrix (its off-diagonal gives the first-order parameter correlations — see correlation). For an unbiased (e.g. linear) estimator the bias is ≈ 0 and the variance ≈ the HESSE error².

jackknife(refit::Function, data; kws...) -> JackknifeResult

Generic (interface-iii) jackknife. data supports length and data[idxvec]; refit(subdata) -> θ̂::AbstractVector. Keywords: d, threaded, and names (parameter labels, default p1, p2, …).

jackknife(cost::AbstractCost, start; kws...) -> JackknifeResult

Delete-1 (or delete-d block) jackknife for a cost object (interface i). Supported for the point-level costs LeastSquares, UnbinnedNLL, ExtendedUnbinnedNLL; binned costs and CostSum raise an ArgumentError. Keywords (d, warm_start, threaded) and the returned JackknifeResult match the model + Data method.

BootstrapResult

Result of bootstrap. Fields:

• samples::Matrix{Float64} — nresample × npar matrix of re-fitted parameter vectors θ̂ (one row per resample, columns in external-parameter order). Rows for non-converged re-fits are kept (as NaN) and flagged in valid.
• valid::Vector{Bool} — per-resample convergence mask (m.valid). Summary statistics below are computed over the valid rows only (see filter_invalid in bootstrap).
• names::Vector{String} — external parameter names.
• estimate::Vector{Float64} — the original full-data optimum θ̂_full (the anchor the re-fits are warm-started from).
• mean::Vector{Float64} — bootstrap mean of θ̂ over valid resamples.
• std::Vector{Float64} — bootstrap standard error (the spread of θ̂; comparable to the HESSE error for a well-specified Gaussian fit).
• ci_lower, ci_upper::Vector{Float64} — percentile confidence interval at ci_level (the [(1-ci_level)/2, (1+ci_level)/2] quantiles of θ̂). These are generally asymmetric about estimate — that asymmetry is the whole point versus the symmetric HESSE error.
• ci_level::Float64 — the CI coverage (default 0.68, i.e. a ±1σ-equivalent).
• covariance::Union{Nothing,Matrix{Float64}} — bootstrap covariance of θ̂ (npar × npar), or nothing unless covariance=true was requested. The off-diagonal captures parameter correlations; correlation returns the standardised correlation matrix from samples regardless of this flag, and the raw samples cloud retains non-Gaussian joint structure a covariance cannot.
• kind::Symbol — :nonparametric (resample points with replacement) or :parametric (regenerate y from the best-fit model + error model).
• nresample::Int, n_valid::Int — total and converged resample counts.
• seed — the master seed used (for reproducibility), or nothing.

JackknifeResult

Result of jackknife. For the delete-1 jackknife there is one re-fit per data point (g = N); for the optional delete-d block jackknife the data are partitioned into g = ceil(N/d) consecutive blocks and one block is deleted per re-fit. Fields:

• samples::Matrix{Float64} — g × npar matrix of leave-one-(group-)out estimates θ̂₍ⱼ₎.
• valid::Vector{Bool} — per-re-fit convergence mask.
• names::Vector{String} — external parameter names.
• estimate::Vector{Float64} — the full-data optimum θ̂_full.
• mean::Vector{Float64} — θ̄ = mean of the leave-one-out estimates.
• bias::Vector{Float64} — jackknife bias estimate (g-1)·(θ̄ - θ̂_full).
• bias_corrected::Vector{Float64} — θ̂_full - bias (the debiased estimate).
• variance::Vector{Float64} — jackknife variance ((g-1)/g)·Σ(θ̂₍ⱼ₎ - θ̄)².
• std::Vector{Float64} — sqrt.(variance) (comparable to the HESSE error).
• covariance::Matrix{Float64} — the full jackknife covariance ((g-1)/g)·Σ(θ̂₍ⱼ₎ - θ̄)(θ̂₍ⱼ₎ - θ̄)ᵀ, whose diagonal is variance. Its off-diagonal captures the parameter correlations (to first order — the jackknife is a linearisation; use bootstrap for non-Gaussian / strongly nonlinear joint structure). See correlation.
• n::Int — number of data points N.
• d::Int — block size (1 for the delete-1 jackknife).
• g::Int — number of deleted groups (N for delete-1).
• n_valid::Int — converged re-fit count.

correlation(r::BootstrapResult) -> Matrix{Float64}
correlation(r::JackknifeResult) -> Matrix{Float64}

The npar × npar Pearson correlation matrix of the parameter estimators, computed from the re-fitted θ̂ over the converged-and-finite resamples (for the jackknife this is the standardised r.covariance). The off-diagonal entries are the estimated parameter correlations — the joint information the per-parameter marginal CIs / errors drop.

Correlation is a linear (second-moment) summary. For a non-Gaussian or nonlinearly-correlated joint distribution (e.g. a curved degeneracy / "banana" in an amplitude fit) inspect the raw joint cloud r.samples directly — a scatter of r.samples[:, i] vs r.samples[:, j], a 2-D density, or a contour — which retains the structure a correlation coefficient cannot. A fixed (zero- variance) parameter yields NaN in its row/column (correlation undefined).

correlation(m::Minuit; skip_fixed=true) -> Matrix{Float64}

The parameter correlation matrix of the fit — the HESSE covariance normalised to unit diagonal (C[i,j] = V[i,j]/√(V[i,i]·V[j,j])). IMinuit.jl/iminuit-compatible (equivalent to matrix(m; correlation=true), and to iminuit's m.covariance.correlation()). Returns nothing if MIGRAD has not run. By default the free-parameter block is returned (skip_fixed=true); pass skip_fixed=false for the full external matrix (fixed parameters then give a 0 row/column).

C = correlation(m)
C[i, j]          # estimated correlation of parameters i and j

find_solution_modes(samples, m::Minuit; whiten=:auto, method=:components,
                    threshold=1.0, min_size=1, min_neighbors=1,
                    fvals=nothing, refine=false, refine_iterate=5,
                    refine_maxfcn=nothing, refine_strategy=nothing,
                    refine_tol=nothing, refine_callback=nothing,
                    parallel=nothing) -> SolutionModes

Cluster a set of Δχ²-accepted parameter samples into DISTINCT solution modes — a "beyond iminuit" capability for multi-modal posteriors. See the file header and docs/src/error_analysis.md for the full rationale.

Arguments

• samples::AbstractMatrix — accepted samples, one parameter vector per row. Width must be either the total parameter count (m.ndim, full external vectors — what the MC-Δχ² sampler get_contours_samples produces) or the free-parameter count (m.npar; fixed parameters are filled from the fit and the clustering ignores them either way).
• m::Minuit — the converged fit the samples were drawn around. Supplies the fallback whitening metric (covariance / errors), the cost function, and — for refine — the full re-fit configuration.

FCN cost (read this for expensive cost functions)

By default (fvals = nothing) the FCN is evaluated at every sample to pick min-χ² representatives and report delta_fval — N full FCN calls before any output (skipped automatically when clustering finds no modes). On expensive FCNs either pass the per-sample χ² you already have (fvals = <vector>; get_contours_samples keeps them), or choose

• fvals = :lazy — evaluate ONLY the K cluster representatives (whitened-space medoids): K FCN calls, full report (χ², delta_fval, χ²-sorted modes), but the representative is the most central member, not the lowest-χ² one.
• fvals = :none — zero FCN calls: representatives are medoids, fval/delta_fval are NaN, and modes are sorted by refined χ² when refine=true, else by population (largest first).

Keyword arguments

• whiten::Symbol=:auto — distance metric (this choice decides whether the tool works at all on its target inputs):
  • :sample — per-coordinate robust cloud scale, σ_k = 1.4826·MAD of the sample column. Fit-independent: it measures the cloud with its own yardstick, which is what a MULTI-BASIN cloud needs — basin separations stay several whitened units wide while intra-basin distances stay ≪ 1. A coordinate with zero cloud spread (constant column) falls back to the fit σ for that coordinate, with a warning. (Note for high-dimensional, sparsely-sampled clouds: with ≳ 8 free parameters and only a few hundred samples, nearest-neighbour distances inside a single cluster approach 1 cloud-σ, so a sparse cluster can fall apart at threshold=1 — the zero-modes diagnostic below reports the cloud's NN scale and the threshold that would reconnect it.)
  • :auto (default) — picks the metric from the cloud/fit width ratio: :sample when the cloud is wider than the fit's local scale in some coordinate (max_k σ_cloud_k/σ_fit_k > 4 — the multi-basin / cross-basin regime where fit-local metrics isolate every point), otherwise the :cov chain below (a single-basin Δχ² cloud sampled at the fit's own scale keeps the statistically tightest metric). Falls back to :errors when neither the cloud nor the covariance provides a usable scale.
  • :cov — full Mahalanobis using the fit's free-parameter covariance (decorrelating + rescaling). The statistically tightest metric when the samples live at the fit's own scale (a single-basin Δχ² cloud from this exact fit). On a cloud spread over many fit-σ (any cross-basin scan, scatter not generated from this fit's errors) every point looks isolated → 0 modes; use :sample there. Falls back to :errors (with a warning) if the covariance is unavailable, not positive-definite, or not trustworthy (!m.accurate — invalid fit or forced-posdef Hessian).
  • :errors — per-parameter fit scale, z_k = x_k / σ_k. No correlations, no matrix inversion; the robust fallback. Same local-metric caveat as :cov. Naive unwhitened Euclidean is intentionally not offered — it is dominated by the largest-scale parameter and merges modes that differ only in a small-scale parameter.
• method::Symbol=:components — :components (built-in single-linkage connected components, no dependency) or :dbscan (density-based; requires using Clustering).
• threshold::Real=1.0 — connection radius in WHITENED units: fit-σ for :cov/:errors, cloud-σ for :sample. Two samples closer than this are linked; single-linkage then chains links into components. Smaller = stricter (won't bridge distinct modes, but may split a sparsely-sampled one); larger risks chaining across a gap. For :dbscan this is the ε radius.
• min_size::Integer=1 — clusters with fewer members are dropped as noise. Default 1 keeps every separated region (the accepted set is already χ²- filtered, so even a few-point region is a candidate solution); raise it to suppress scatter.
• min_neighbors::Integer=1 — :dbscan core-point density (ignored by :components).
• fvals::Union{Nothing,Symbol,AbstractVector}=nothing — per-sample χ² policy: a vector of precomputed values, nothing (evaluate the FCN at all samples), :lazy (K representative evaluations) or :none (no FCN calls). See FCN cost above.
• refine::Bool=false — re-run MIGRAD from each cluster's representative to recover that mode's true local minimum + its own errors. Flags new_min if a mode is DEEPER than the global best (the main fit missed the better basin).
• refine_iterate::Integer=5 — iterate passed to each re-fit migrad! (max MIGRAD attempts per mode).
• refine_maxfcn::Union{Nothing,Integer}=nothing — FCN call cap per MIGRAD attempt of a re-fit. Reaching the cap also stops the retry loop, so a budgeted mode costs ≈ 1.3 × refine_maxfcn calls in practice (an attempt may extend its own budget by ×1.3 — the C++ second-pass bump — and the cap is checked once per MIGRAD iteration, so it can overshoot by one iteration's evaluations); refine_maxfcn × refine_iterate is a hard upper bound. nothing = Minuit's default budget 200 + 100·n + 5·n² per attempt — effectively uncapped; SET THIS on expensive FCNs.
• refine_strategy::Union{Nothing,Strategy,Integer}=nothing — MIGRAD strategy for the re-fits (nothing inherits the parent fit's; 0 = fastest, for triage on slow FCNs).
• refine_tol::Union{Nothing,Real}=nothing — EDM tolerance for the re-fits (nothing inherits the parent fit's).
• refine_callback=nothing — a function called once per finished re-fit (in completion order), for checkpointing long runs on slow FCNs: a killed job then loses at most the mode in flight. Receives a NamedTuple (k, K, representative, n_points, member_indices, refined, refined_values, refined_errors, refined_fval, refined_valid, refined_nfcn, new_min, walltime) where k is the completion count (1…K). May be invoked from worker threads when parallel, but invocations are serialized (a lock), so writing to a file/IO from it is safe; exceptions it throws are caught and warned, never aborting the run.
• parallel::Union{Bool,Nothing}=nothing — parallelize FCN evaluation and per-mode re-fits across threads. nothing ⇒ auto-on when m.threaded_gradient is set and Threads.nthreads() > 1 (the same FCN thread-safety contract as Phase G/H threaded gradients). Set true/false to force.

Returns

A SolutionModes (an AbstractVector{SolutionMode}), sorted with the main solution first (lowest χ²; see fvals = :none above for the no-χ² ordering). Pretty-prints a summary report.

When clustering finds no modes (or bins more than half the samples as noise), a diagnostic warning reports the cloud's median nearest-neighbour whitened distance against threshold and suggests the concrete fix (a cloud-scaled metric or a larger threshold) — a 0-mode result on a sane cloud almost always means the metric, not the cloud.

Example

m = Minuit(chi2, x0; names=pnames); migrad!(m)
samples = get_contours_samples(m; ...)          # MC-Δχ² accepted set (rows = vectors)
modes = find_solution_modes(samples, m)          # cluster into distinct solutions
length(modes) > 1 && @warn "multi-modal: $(length(modes)) distinct solutions"
modes_refined = find_solution_modes(samples, m; refine=true)   # + per-mode re-fit

# Expensive FCN (seconds per call): zero-eval triage, then budgeted refine
# with per-mode checkpointing.
modes = find_solution_modes(samples, m; fvals = :none)
modes = find_solution_modes(samples, m; fvals = :lazy, refine = true,
                            refine_maxfcn = 500, refine_strategy = 0,
                            refine_callback = r -> println("mode $(r.k)/$(r.K): ",
                                                           r.refined_fval))

SolutionMode

One distinct solution found by find_solution_modes: a connected cluster of Δχ²-accepted samples that is separated (in whitened parameter space) from the other clusters.

Fields

• index::Int — rank, 1 = main solution (lowest χ²), 2, 3, … by ascending χ². (With fvals = :none the per-sample χ² is unknown: modes are ranked by refined χ² when refine=true, else by cluster population, largest first.)
• representative::Vector{Float64} — the minimum-χ² sample in the cluster (full external parameter vector, length = total parameter count). With fvals = :none or :lazy the representative is instead the cluster's whitened-space medoid (most central member).
• fval::Float64 — χ² (cost) at representative (NaN with fvals = :none).
• delta_fval::Float64 — fval minus the global-best χ². ≥ 0 for every mode when the global best is the fit minimum; the main mode is ≈ 0. Interpret against the fit's errordef (up): a mode with delta_fval ≲ up is statistically comparable to the main solution. (NaN with fvals = :none.)
• param_ranges::Vector{Tuple{Float64,Float64}} — per-parameter (min, max) over the cluster's samples (full external coordinates).
• n_points::Int — number of samples in the cluster.
• fraction::Float64 — n_points / total accepted samples.
• member_indices::Vector{Int} — row indices (into the samples matrix) of this cluster's members.

Re-fit fields, populated only when find_solution_modes(...; refine=true):

• refined::Bool — whether a per-mode MIGRAD re-fit was run and succeeded.
• refined_values::Vector{Float64} — re-fit parameter values (empty if not refined).
• refined_errors::Vector{Float64} — re-fit errors (empty if not refined).
• refined_fval::Float64 — re-fit χ² (NaN if not refined).
• refined_valid::Bool — whether the re-fit MIGRAD validated.
• refined_nfcn::Int — FCN calls used by the re-fit.
• new_min::Bool — true if this mode's re-fit reached a DEEPER minimum than the global best — i.e. the main fit missed the better basin. Flagged prominently in the report; connects to the IAM cold-start convergence gap (a separated cluster can be the basin the global fit failed to find).
• refined_walltime::Float64 — wall-clock seconds the re-fit took (NaN if refine=false).

SolutionModes <: AbstractVector{SolutionMode}

Result of find_solution_modes: an indexable, iterable vector of SolutionModes (so modes[1], length(modes), for s in modes all work) plus the metadata needed for the summary report. Pretty-prints a "Found K distinct solutions within Δχ²" table; index it to get the individual modes.

Fields (besides the modes themselves)

• global_best::Float64 — reference χ² the delta_fvals are measured from.
• up::Float64 — the fit's errordef (1.0 χ², 0.5 NLL) — the natural Δχ² yardstick.
• whiten::Symbol — the metric actually used (:sample, :cov or :errors; may differ from the request after a fallback, and is never :auto — :auto resolves to :sample or :errors).
• method::Symbol — :components or :dbscan.
• threshold::Float64 — whitened connection radius (σ units).
• n_noise::Int — samples dropped as noise (clusters smaller than min_size).
• n_samples::Int — total accepted samples clustered.
• param_names::Vector{String} — parameter names (for the report).

find_deeper_minimum(m::Minuit; kwargs...) -> Minuit

Parameter-perturbation basin-hopping search for a deeper minimum, starting from the fit m. Each round draws n_restarts perturbed restarts around the current best (each FREE coordinate jittered by perturb · scaleᵢ · randn, with scaleᵢ = max(|xᵢ|, |errorᵢ|, abs_floor)), MIGRADs each, and adopts any deeper valid minimum; it stops when a round finds no improvement or after max_rounds.

Every restart is a clone of m (via Minuit), so m's parameter limits and fixed parameters are honoured — fixed parameters are never jittered and stay pinned at their value, and jittered coordinates are clamped into their bounds. The search therefore stays inside the same constrained parameter space as your fit. Returns the deepest Minuit found (MIGRAD + HESSE already run); check its .valid property. m itself is not mutated.

Keyword arguments

• n_restarts::Integer = 24 — perturbed restarts per round (≥ 1).
• perturb::Real = 1.0 — exploration radius, as a multiple of each parameter's scale. The key knob on a hard surface.
• abs_floor::Real = 0.0 — absolute lower bound on each coordinate's jitter scale (raise it for a parameter sitting near 0 with a tiny step).
• max_rounds::Integer = 50 — safety backstop, not the normal stop. The search stops when a round finds no deeper basin (convergence); max_rounds only bounds a pathological non-converging run. A @warn fires if it is hit while still improving.
• strategy = m.strategy — MIGRAD strategy for every (re-)fit.
• maxfcn::Union{Integer,Nothing} = nothing — per-fit MIGRAD call budget.
• min_improvement::Real = 1e-3 — χ² drop required to adopt a restart.
• seed::Union{Integer,Nothing} = nothing — RNG seed for reproducible restarts.
• verbose::Bool = false — log the best χ² per round.

See also find_solution_modes and the data-resampling overload find_deeper_minimum(m, refit, data).

find_deeper_minimum(cf, x0, errors; limits=nothing, fixed=nothing, names=nothing, kwargs...) -> Minuit

Parameter-perturbation search from a cost function / callable cf and start x0/errors. Builds a Minuit with the given limits, fixed flags and names (same meaning as the Minuit constructor) and delegates to find_deeper_minimum(m::Minuit) — so bounds and fixed parameters are honoured. cf may be an AbstractCostFunction (its gradient, if any, is carried) or a bare callable (up sets the error definition).

find_deeper_minimum(m::Minuit, refit, data; kwargs...) -> Minuit

Data-resampling basin-hopping search for a deeper minimum on a multi-basin data-fitting objective. Each round bootstrap-resamples data and re-fits each resample with refit (those drift toward whichever basin best explains that subset); the candidates are clustered with find_solution_modes(...; whiten=:cov, refine=true) (the anchor-fit metric — converged candidates are tight clumps at the fit's own scale) and re-evaluated on the ORIGINAL objective (so the χ² comparison is honest); the deepest valid new basin is adopted, and the loop repeats.

Constraints are honoured. The re-fit/refinement runs through Minuit(m.fcn.f, m) (which carries m's limits + fixed flags and pins fixed parameters), and the adoption rebuild clones the same configuration — so a fixed parameter stays fixed and a bounded one stays in bounds throughout. Returns the deepest Minuit (MIGRAD

• HESSE already run); check .valid. m is not mutated.

Arguments

• m::Minuit — a converged fit. Supplies the cost function, current best as the warm start, the HESSE covariance for whitening, and the constraint structure.
• refit — refit(subdata, start::Vector{Float64}) -> Vector{Float64} (a NaN- filled vector of the same length ⇒ invalid, dropped). Accepts functors. To keep the discovery itself constrained, refit should apply the same fix/limits; any fixed parameter is in any case re-pinned during refinement.
• data — the full dataset, supporting data[idx] bootstrap indexing.

Keyword arguments

• n_discovery::Integer = 20 — bootstrap resamples per round (≥ 2; keep ≥ 10 for stable clustering).
• max_rounds::Integer = 50 — safety backstop, not the normal stop. The search stops when a round adopts no deeper basin (convergence); a @warn fires if the cap is hit while still improving.
• strategy = m.strategy — MIGRAD strategy for the adoption re-fit.
• min_improvement::Real = 1e-3 — χ² drop required to adopt a basin.
• parallel::Union{Bool,Nothing} = nothing — true threads the discovery loop (only if refit is thread-safe) AND is passed to find_solution_modes; nothing/false runs discovery serially.
• seed::Union{Integer,Nothing} = nothing — RNG seed (indices are pre-generated sequentially, so results are deterministic regardless of parallel).
• verbose::Bool = false — log χ² improvement per round.

Suitability check

If the first round finds no deeper basin, a @warn is emitted and a fitted clone of m (the same minimum, since m is never mutated) is returned — the surface may be single-basin here, or bootstrap coverage was insufficient (raise n_discovery); for parameter-space search try the perturbation overload find_deeper_minimum(m; perturb=…).

find_deeper_minimum(cf, x0, errors, refit, data; limits=nothing, fixed=nothing, names=nothing, kwargs...) -> Minuit

Data-resampling search from a cost function / callable. Builds a Minuit with the given limits/fixed/names, fits it, and delegates to the (m::Minuit, refit, data) overload (constraints honoured). When you already hold a converged Minuit, prefer passing it directly.

optim(m::Minuit; method=:lbfgs, ncall=nothing, maxcall=nothing,
                 tol=nothing, options=nothing) -> Minuit

The Julia analog of iminuit's m.scipy(...) — the alternative-minimizer escape hatch, powered by Optim.jl instead of Python's scipy.optimize. Load using Optim to enable (it is an optional package extension).

Minimises the FCN with the chosen Optim optimizer starting from m's current parameter values, honours fixed parameters and box limits (via Optim's Fminbox), writes the optimum back into m (so m.values / m.fval update), and returns m. Run hesse(m) afterwards for the covariance — matching iminuit's scipy-then-hesse flow. Use this when MIGRAD struggles (stiff / ill-conditioned problems where a trust-region or derivative-free method does better).

method mapping (case / dash / underscore insensitive)

Derivative-free (:neldermead) and second-order (:newton) methods cannot be combined with box limits — Optim's Fminbox requires a first-order optimizer. Use a first-order method (:lbfgs / :bfgs / :conjugategradient / :gradientdescent) for bounded fits. When the constructor was given grad=…, the analytical gradient is passed through to first-order optimizers; otherwise Optim finite-differences it.

Keyword arguments

• method — optimizer selector (see table). Default :lbfgs.
• ncall / maxcall — function-evaluation budget → Optim's f_calls_limit.
• tol — gradient-norm convergence tolerance → Optim's g_tol.
• options — a full Optim.Options(...) for fine control (overrides ncall/maxcall/tol).

For full control over the optimizer object itself, use minimize_with: minimize_with(m, Optim.LBFGS()).

Examples

using NativeMinuit, Optim
m = Minuit(fcn, x0)
optim(m; method=:lbfgs)   # or m |> optim
hesse(m)                  # covariance, à la iminuit

minimize_with(m::Minuit, optimizer=nothing; method=:lbfgs, ncall=nothing,
              maxcall=nothing, tol=nothing, options=nothing) -> Minuit

Clearer-named alias of optim: minimise m's FCN with an alternative optimizer from Optim.jl, write the optimum back, and return m (run hesse(m) afterwards for the covariance). Load using Optim to enable.

The first positional argument may be an Optim optimizer object, bypassing the method name table for full control:

using NativeMinuit, Optim
minimize_with(m, LBFGS())                       # optimizer object
minimize_with(m, NelderMead())
minimize_with(m; method=:bfgs, tol=1e-10)       # by name, like optim(m; …)

See optim for the method mapping and keyword semantics.

to_latex(m::Minuit; siunitx=true, booktabs=true,
         caption=nothing, label=nothing) -> String

Render the fitted parameters of m as a publication-ready LaTeX table.

Defaults to a booktabs rule set with siunitx \num{} numbers. Asymmetric MINOS errors (when minos! has run) are written as \num{x}^{+hi}_{-lo}; otherwise a symmetric \num{x} \pm \num{e} is used. Numbers are rounded to the uncertainty (1–2 significant figures on the error). Fixed parameters show the value alone, tagged (fixed).

Set siunitx=false for plain numbers (no \num{}), booktabs=false for \hline rules, and pass caption/label to wrap the tabular in a table float.

Requires \usepackage{booktabs} and \usepackage{siunitx} in the document preamble (unless the corresponding option is disabled).

Example

m = Minuit(p -> (p[1]-1)^2 + (p[2]-2)^2, [0.0, 0.0]; names=["mass", "width"])
migrad!(m)
print(to_latex(m))

emits (numbers rounded to the uncertainty; trailing zeros stripped):

\begin{tabular}{l c}
\toprule
Parameter & Value \\
\midrule
mass & $\num{1} \pm \num{1}$ \\
width & $\num{2} \pm \num{1}$ \\
\bottomrule
\end{tabular}

to_latex(e::MinosError; value=e.min_par_value, siunitx=true) -> String

Render a single asymmetric MINOS result as \num{value}^{+hi}_{-lo} (no surrounding math delimiters), suitable for dropping into running text.

If either MINOS side failed to converge (a non-finite upper/lower), there is no meaningful asymmetric error, so the central value alone is returned (\num{value}).

mn_plot_text(pts::AbstractVector{<:Tuple{Real,Real}};
             width=60, height=20,
             par_x="x", par_y="y",
             x_center=nothing) -> String

Lower-level ASCII renderer for an arbitrary 2D point set (e.g. the raw Vector{Tuple{Float64,Float64}} returned by mncontour). Returns a multi-line String ready for println / @info.

The bounding box auto-scales to the data with a small margin and snaps to round-number ticks via the Minuit2 mnbins heuristic. Cells are single characters: * for a single point, & where two differing characters collide (matches mntplot.cxx semantics — same-character stamps coalesce silently). x_center, when supplied, marks the minimum with X.

Width / height are hints — the actual grid size is the nearest "nice" binning, so the rendered box may be a little larger than requested.

Defensive behavior

• Empty pts → returns "EMPTY plot — no points to render.\n".
• Non-finite (NaN / ±Inf) points are silently skipped; if no finite points remain, behaves like the empty case.

mn_plot_text(c::ContoursError; width=60, height=20) -> String

ASCII-render a 2-parameter MnContours result as a Cartesian box. Returns a single multi-line String suitable for println. The fitted minimum (from the embedded MINOS errors) is marked with X.

Invalid contours (c.valid == false or no points) render as an empty-plot message rather than throwing.

draw_contour(m::Minuit, par1, par2; size=50, bound=2, kws...) -> Plots.Plot

Filled-contour plot of the FCN grid slice from contour_grid(m, par1, par2; subtract_min=true) — iminuit's m.draw_contour. A landscape view, NOT a confidence region (see contour_grid); for confidence contours use draw_mncontour. Requires using Plots. (≤ 0.4 this drew the contour_ellipse approximation instead.)

draw_mncontour(m::Minuit, par1, par2; numpoints=100, cl=nothing, kws...) -> Plots.Plot

Plot exact MINOS 2D confidence contours from mncontour (boundary search with per-point re-minimization) at one or several confidence levels — cl follows mncontour's iminuit semantics (default → joint 2-D 68 % region) and may be a vector to overlay several contours. Requires using Plots. (≤ 0.4 this drew the fast contour_ellipse approximation at Δχ² = up instead.)

draw_profile(m::Minuit, par; bins=100, kws...) -> Plots.Plot

Plot the 1D scan from profile(m, par; bins=bins). Requires using Plots.

draw_mnprofile(m::Minuit, par; bins=30, kws...) -> Plots.Plot

Plot the 1D MINOS profile from mnprofile(m, par; bins=bins). Requires using Plots.

draw_mnmatrix(m::Minuit; kws...) -> Plots.Plot

Plot all pairwise 2D contours arranged as a triangular matrix. Requires using Plots.