def calculate_peaking_factor(central: float,
                           volume_avg: float) -> float:
    r"""Peaking factor, central value over volume average.

    $$\mathrm{PF} = \frac{X(0)}{\langle X\rangle}$$

    Parameters
    ----------
    central : float
        Value on the magnetic axis [any].
    volume_avg : float
        Volume average of the same quantity, same unit [any].

    Returns
    -------
    float
        Peaking factor [-].

    Numerical notes
    ---------------
    A zero volume average warns and returns ``nan``.

    Reduction
    ---------
    input: scalar_0d
    output: scalar_0d
    kind: dimensionless_normalization
    locality: global
    role: profile_descriptor
    """
    return _guarded_ratio(
        central, volume_avg, what="calculate_peaking_factor", because="volume_avg"
    )

def calculate_poloidal_flux(R: np.ndarray,
                          B_theta: np.ndarray,
                          l: np.ndarray,
                          psi_axis: float = 0.0) -> float:
    r"""Poloidal flux per radian from a line integral of $RB_\theta$.

    $$\psi(l) = \int_0^{l} R\,B_\theta\,dl' + \psi_a$$

    Parameters
    ----------
    R : np.ndarray
        Major radius along the path [m].
    B_theta : np.ndarray
        Poloidal field normal to the path [T].
    l : np.ndarray
        Path coordinate, monotonic [m].
    psi_axis : float, optional
        Offset added to the integral; default 0 [Wb/rad].

    Returns
    -------
    float
        Flux at the end of the path [Wb/rad].

    Convention
    ----------
    Flux per radian (COCOS 1-8 storage); multiply by $2\pi$ for the IMAS
    full-weber flux.  Sign follows ``B_theta`` and the direction of ``l``.
    The physics wrapper is :func:`vaft.formula.equilibrium.psi_from_RBtheta`.

    Numerical notes
    ---------------
    Trapezoidal rule over the whole path; returns the end value only.
    """
    return trapz_integral(l, R * B_theta) + psi_axis

def calculate_toroidal_flux(B_phi: np.ndarray,
                            R: np.ndarray = None,
                            Z: np.ndarray = None,
                            mask: np.ndarray = None,
                            *,
                            dA: np.ndarray = None) -> float:
    r"""Toroidal flux $\Phi = \int B_\varphi\,dR\,dZ$ over a rectilinear $(R, Z)$ grid.

    $$\Phi = \iint_{S} B_\varphi(R, Z)\,dR\,dZ$$

    by the 2-D trapezoidal rule, over the whole grid or over the cells a
    ``mask`` selects.

    Parameters
    ----------
    B_phi : np.ndarray
        Toroidal field on the grid, indexed ``[R, Z]`` (IMAS ``profiles_2d``
        order) [T].
    R : np.ndarray, optional
        Major radius: a 1-D strictly monotonic axis of length ``B_phi.shape[0]``,
        or an ``ij``-indexed 2-D mesh of the shape of ``B_phi`` [m].
    Z : np.ndarray, optional
        Height: a 1-D strictly monotonic axis of length ``B_phi.shape[1]``, or an
        ``ij``-indexed 2-D mesh of the shape of ``B_phi`` [m].
    mask : np.ndarray, optional
        Boolean map of the shape of ``B_phi``; ``True`` marks the samples inside
        the cross-section.  ``psi_N <= 1`` alone also selects the private-flux
        region below a diverted X-point; combine it with an inside-the-LCFS
        test there.  Default: the whole grid [-].
    dA : np.ndarray, optional
        Deprecated (#358).  Caller-supplied area elements of the old bare sum
        $\sum_i B_{\varphi,i}\,\Delta A_i$; also taken from a two-argument call
        ``calculate_toroidal_flux(B_phi, dA)`` when ``B_phi`` is not 2-D or ``dA``
        is a scalar cell area (a 2-D area map must pass ``dA=`` by keyword).  Emits a ``DeprecationWarning`` [m^2].

    Returns
    -------
    float
        Toroidal flux through the selected area [Wb].

    Raises
    ------
    ValueError
        If the grid is not rectilinear, an axis is not strictly monotonic, a
        shape does not match ``B_phi``, or ``mask`` is not boolean.
    TypeError
        If neither the grid nor ``dA`` is given, both are, or a 2-D map is
        given a second argument without ``Z``.

    Convention
    ----------
    Full weber (toroidal flux has no per-radian form); sign of $B_\varphi$
    along the COCOS $\varphi$ direction, so a reversed-field discharge flips it
    as surely as a COCOS conversion does.  The area is counted positive
    whichever direction the axes are stored in, so reversing ``R`` or ``Z``
    does not flip the sign.  The physics wrapper is
    :func:`vaft.formula.equilibrium.phi_from_Bphi`.

    Numerical notes
    ---------------
    Trapezoidal rule along ``Z`` then ``R`` through
    :func:`vaft.compat.trapz_compat`: second order in the spacing over the full
    grid, exact for a field bilinear in $(R, Z)$.  With a ``mask`` the integrand
    is set to zero outside it, so the staircase boundary makes the error first
    order in the cell size; refine the grid, not the mask, to converge.  The
    deprecated ``dA`` path is the old Riemann sum, first order and only as good
    as the caller's areas.

    References
    ----------
    .. [1] J. Wesson, *Tokamaks*, 4th ed., Oxford University Press (2011),
           Sec. 3.4 (toroidal flux and the safety factor).
    """
    return _toroidal_flux(B_phi, R, Z, mask, dA, stacklevel=3)

def calculate_volume_weighted_average(x: np.ndarray,
                                    V: np.ndarray) -> float:
    r"""Volume-weighted average of a sampled profile.

    $$\langle X\rangle = \frac{\sum_i X_i\,V_i}{\sum_i V_i}$$

    Parameters
    ----------
    x : np.ndarray
        Profile values at each cell [any].
    V : np.ndarray
        Volume of each cell, same shape [m^3].

    Returns
    -------
    float
        Volume-weighted average in the unit of ``x`` [any].

    Numerical notes
    ---------------
    A zero total volume warns and returns ``nan``.
    """
    return _guarded_ratio(
        np.sum(x * V),
        np.sum(V),
        what="calculate_volume_weighted_average",
        because="sum(V)",
    )

def eped_tanh_bounds(x, y):
    """Parameter box the ``eped_tanh`` fit uses for this data.

    Public so that a caller checking whether a fitted parameter came back
    resting on a bound compares against the same box the fit used, rather than
    re-deriving one that can drift out of step with it.  The amplitude bounds
    scale with the data's magnitude; the shape bounds are fixed.

    Parameters
    ----------
    x : array_like
        Radial coordinate; used only for its length [-].
    y : array_like
        Profile samples, whose magnitude sets the amplitude bounds [any].

    Returns
    -------
    tuple of list
        ``(lower, upper)``, each seven entries ordered as the model's
        arguments ``f0, f_ped, f_sep, x_ped, width, alpha, beta``.  The
        amplitudes are ``+-5 max|y|`` [any]; ``x_ped`` is in ``[0.75, 0.999]``
        and ``width`` in ``[0.01, 0.3]`` [-]; ``alpha`` and ``beta`` are in
        ``[1, 8]`` [-].

    Notes
    -----
    ``alpha`` and ``beta`` stay at or above one because below it the gated
    core term has an infinite slope -- at ``x_ped`` for ``beta``, at the
    magnetic axis for ``alpha`` -- which no flux-function profile has, and
    which would make any gradient taken from the fit depend on the caller's
    grid spacing.

    References
    ----------
    .. [EPED] The seven-parameter form and the ``x_ped``/``width`` box follow
       the pedestal-fitting study ported in migration decision D-05;
       :func:`vaft.process.profile.pedestal_top` is the consumer.
    """
    _, lower, upper = _initial_eped_tanh_guess(x, y)
    return lower, upper

def gp_fit(x, y, y_std, x_eval, *, n_restarts_optimizer: int = 5, random_state: int = 0):
    """Fit a squared-exponential Gaussian process and evaluate it on a grid.

    Parameters
    ----------
    x : array-like
        Data abscissae, 1-D [any].
    y : array-like
        Data values [any].
    y_std : array-like or None
        Per-point one-sigma uncertainty, entering as the diagonal noise;
        ``None`` uses a numerical jitter [any].
    x_eval : array-like
        Evaluation grid [any].
    n_restarts_optimizer : int, optional
        Extra random restarts of the marginal-likelihood optimisation, on top of
        the one from the default hyperparameters; default 5 [-].
    random_state : int, optional
        Seed for the restart draws, so a fit is reproducible; default 0 [-].

    Returns
    -------
    mean : np.ndarray
        Posterior mean on ``x_eval`` [any].
    std : np.ndarray
        Posterior standard deviation on ``x_eval`` [any].

    Numerical notes
    ---------------
    ``y`` is standardised before fitting and the posterior is mapped back,
    which is what makes one set of hyperparameter bounds work for data in
    m^-3 and in keV alike.  The marginal likelihood is optimised with
    ``scipy.optimize.minimize`` (L-BFGS-B) over the log-hyperparameters.

    References
    ----------
    .. [1] C. E. Rasmussen and C. K. I. Williams, *Gaussian Processes for
           Machine Learning*, MIT Press (2006), Algorithm 2.1 and Eq. (5.8).
    """
    state = _gp_train(
        x, y, y_std,
        n_restarts_optimizer=n_restarts_optimizer,
        random_state=random_state,
    )
    return _gp_evaluate(state, x_eval)

def gradient(x: np.ndarray, y: np.ndarray) -> np.ndarray:
    r"""Derivative $dy/dx$ along the first axis of a sampled profile or map.

    $$\frac{dy}{dx}\Big|_i \approx \frac{y_{i+1} - y_{i-1}}{x_{i+1} - x_{i-1}}$$

    (second-order central difference on a possibly non-uniform grid).

    Parameters
    ----------
    x : np.ndarray
        Independent variable, monotonic [any].
    y : np.ndarray
        Dependent variable whose first axis has the length of ``x``: a 1-D
        profile or an N-D map indexed ``[x, ...]`` [any].

    Returns
    -------
    np.ndarray
        $dy/dx$ at every sample, of the shape of ``y``, in units of ``y`` per
        unit ``x`` [any].

    Numerical notes
    ---------------
    Wraps ``numpy.gradient(y, x, axis=0)``: second-order accurate in the
    interior, first-order one-sided at the two end points, and
    noise-amplifying (a relative noise $\delta$ on ``y`` becomes
    $\delta/\Delta x$ on the derivative).  Needs at least two samples; only
    the first axis of an N-D ``y`` is differentiated, each column
    independently, so the result equals stacking the 1-D derivatives.
    """
    return np.gradient(y, x, axis=0)

def make_fit_function(mode):
    r"""Build a 1-D parametric model $f(x; c_0, c_1, \dots)$ for profile fitting.

    $$\begin{aligned}
      \text{polynomial:}\ & (1-x)\,\textstyle\sum_k c_kx^k &
      \text{free\_polynomial:}\ & \textstyle\sum_k c_kx^k \\
      \text{exponential:}\ & (1-x)\exp\big(\textstyle\sum_k c_kx^k\big) &
      \text{free\_exponential:}\ & \exp\big(\textstyle\sum_k c_kx^k\big)
    \end{aligned}$$

    Parameters
    ----------
    mode : str
        Model name, case-insensitive [str].
        One of ``'polynomial'``, ``'free_polynomial'``, ``'exponential'``,
        ``'free_exponential'`` (a few aliases are accepted).

    Returns
    -------
    callable
        ``f(x, *coeffs)`` evaluating the model, in the unit of the data [any].

    Raises
    ------
    ValueError
        For an unknown mode.

    Assumptions
    -----------
    ``x`` is a normalised radius on $[0, 1]$: the $(1 - x)$ factor forces the
    constrained modes to zero at $x = 1$; the exponential modes are strictly
    positive and decay monotonically when the polynomial is decreasing.
    """
    mode = mode.lower()
    if mode == 'polynomial':
        # (1-x)*poly(x) -> enforces value -> 0 at x=1
        def func(x, *coeffs):
            x = np.asarray(x, dtype=float)
            s = 0.0
            for k in range(len(coeffs)):
                s = s + coeffs[k] * x**k
            return (1.0 - x) * s
    elif mode in {'free_polynomial', 'polynomial_unconstrained', 'unconstrained_polynomial'}:
        # plain poly(x) -> no boundary constraint at x=1
        def func(x, *coeffs):
            x = np.asarray(x, dtype=float)
            s = 0.0
            for k in range(len(coeffs)):
                s = s + coeffs[k] * x**k
            return s
    elif mode == 'exponential':
        # (1-x)*exp(poly(x)) -> goes to 0 at x=1, stays positive
        def func(x, *coeffs):
            x = np.asarray(x, dtype=float)
            s = 0.0
            for k in range(len(coeffs)):
                s = s + coeffs[k] * x**k
            return (1.0 - x) * np.exp(s)
    elif mode in {'free_exponential', 'exp_free', 'exponential_unconstrained'}:
        # exp(poly(x)) -> always > 0, NO edge-zero; monotonic decay if poly decreasing.
        # Edge value exp(poly(1)) stays small-but-finite (never exactly 0, never rises if c1<0).
        def func(x, *coeffs):
            x = np.asarray(x, dtype=float)
            s = 0.0
            for k in range(len(coeffs)):
                s = s + coeffs[k] * x**k
            return np.exp(s)
    else:
        raise ValueError(f"Invalid fitting function: {mode}")
    return func

def normalize_profile(x: Union[float, np.ndarray],
                     x_axis: float,
                     x_boundary: float) -> Union[float, np.ndarray]:
    r"""Linear normalisation of a profile between its axis and boundary values.

    $$x_N = \frac{x - x_{\mathrm{axis}}}{x_{\mathrm{boundary}} - x_{\mathrm{axis}}}$$

    Parameters
    ----------
    x : float or np.ndarray
        Values to normalise [any].
    x_axis : float
        Value mapped to 0 [any].
    x_boundary : float
        Value mapped to 1 [any].

    Returns
    -------
    float or np.ndarray
        Normalised values, 0 at the axis value and 1 at the boundary value [-].

    Numerical notes
    ---------------
    A degenerate profile whose axis and boundary values are equal -- which
    packaged VEST samples do contain -- warns and returns ``nan`` rather than
    propagating ``inf``.
    """
    return _guarded_ratio(
        x - x_axis,
        x_boundary - x_axis,
        what="normalize_profile",
        because="x_boundary - x_axis",
    )

def normalized_gradient_scale_length(x: np.ndarray, y: np.ndarray, a: float) -> np.ndarray:
    r"""Normalised inverse gradient scale length $a/L_y$ of a sampled profile.

    $$\frac{a}{L_y} = -\frac{a}{y}\,\frac{dy}{dx}$$

    Parameters
    ----------
    x : np.ndarray
        Radial coordinate, 1-D and monotonic, e.g. minor radius $r$ [m].
    y : np.ndarray
        Profile sampled on ``x`` along its last axis, e.g. $T_e$ or $n_e$;
        shape ``(..., len(x))`` [any].
    a : float
        Normalising length, in the unit of ``x``; for a gyrokinetic input the
        minor radius [m].

    Returns
    -------
    np.ndarray
        $a/L_y$ at every sample, of the shape of ``y``; ``nan`` where ``y``
        is zero or non-finite, with a ``RuntimeWarning`` [-].

    Convention
    ----------
    Positive for a profile that falls with ``x`` -- the usual peaked core
    profile -- so a larger number is a steeper drive.  ``x`` and ``a`` must
    share a unit and ``x`` must be the radius $a$ normalises: on $\rho_N$ or
    $\sqrt{\psi_N}$ pass ``a = 1`` and the result is $-d\ln y/d\rho$, which
    differs from $a/L_y$ by the Jacobian $d\rho/dr$.  The same holds for
    $a/L_{T_e}$, $a/L_{n_e}$ and $a/L_{p}$: only the profile changes.

    Numerical notes
    ---------------
    The derivative is ``numpy.gradient(y, x, axis=-1)``, identical to
    :func:`gradient` for a 1-D profile: second-order central differences in
    the interior, first-order one-sided at the two ends.  Differentiation
    amplifies noise, and dividing by ``y`` amplifies it again where the
    profile is small -- the edge -- so smooth or fit the profile first.

    Reduction
    ---------
    input: profile_1d
    output: profile_1d
    kind: differential
    locality: flux_surface_local
    role: closure_input

    References
    ----------
    .. [1] G. M. Staebler, J. E. Kinsey and R. E. Waltz, Phys. Plasmas 14
           (2007) 055909 (the $a/L_T$, $a/L_n$ drive inputs of TGLF).

    See Also
    --------
    gradient
    """
    x = np.asarray(x, dtype=float)
    y = np.asarray(y, dtype=float)
    derivative = np.gradient(y, x, axis=-1)
    return _guarded_ratio(
        -float(a) * derivative,
        y,
        what="normalized_gradient_scale_length",
        because="the profile value y",
    )

def trapz_integral(x: np.ndarray, y: np.ndarray) -> float:
    r"""Definite integral $\int y\,dx$ by the trapezoidal rule.

    $$\int y\,dx \approx \sum_i \frac{y_i + y_{i+1}}{2}\,(x_{i+1} - x_i)$$

    Parameters
    ----------
    x : np.ndarray
        Sample abscissae, monotonic [any].
    y : np.ndarray
        Integrand at the samples [any].

    Returns
    -------
    float
        Integral over the sampled range, in units of ``y`` times ``x`` [any].

    Numerical notes
    ---------------
    ``numpy.trapezoid`` through :func:`vaft.compat.trapz_compat`; second-order
    accurate in the spacing, exact for piecewise-linear integrands, and
    sign-reversed for a decreasing ``x``.  A same-named helper in
    :mod:`vaft.formula.green` shadows this one on the package namespace
    (``vaft.formula.trapz_integral`` is the Green's-function copy).
    """
    return trapz_compat(y, x=x)