def evaluate_regularized_profile(psi_n, parameters, config, *, derivative=0,
                                 allow_extrapolation=False):
    r"""Evaluate a positive regular flux profile or its flux derivative.

    $$F(x)=s\sum_j B_j(x/d)\,[f+\log(1+\exp(\theta_j))]$$

    Parameter order and basis [1]_ are shared with the offline fit. Finite
    F' implies dF(rho_pol squared)/drho_pol=0 on axis without fixing F'(0).

    Parameters
    ----------
    psi_n : array_like
        Finite normalized poloidal flux coordinates [-].
    parameters : array_like
        Finite unconstrained softplus parameters in basis order [-].
    config : RegularizedProfileConfig
        Basis, floor and domain declaration, with physical scale s [-].
    derivative : int, optional
        Flux derivative order, zero, one or two [-].
    allow_extrapolation : bool, optional
        Enable diagnostic evaluations outside the positivity domain [-].

    Returns
    -------
    values : ndarray
        Profile or derivative in the physical unit y of config.scale [y].

    Limitations
    -----------
    Outside the declared domain the default is NaN. Explicit extrapolation
    carries no positivity guarantee. Numerical overflow is not clipped.

    References
    ----------
    .. [1] C. de Boor, A Practical Guide to Splines, revised edition (2001),
       Chapters IX-X, B-spline representation and differentiation.
    """
    x = np.asarray(psi_n, float)
    theta = np.asarray(parameters, float)
    if theta.shape != (config.coefficient_count,) or not np.isfinite(theta).all():
        raise ValueError("finite parameter vector of the declared size required")
    out = config.scale * (profile_basis(x, config, derivative) @
                          (config.floor + np.logaddexp(0., theta)))
    if not allow_extrapolation:
        out = np.where((x >= 0) & (x <= config.domain_max), out, np.nan)
    return out

def profile_basis(psi_n, config, derivative=0):
    r"""Evaluate the positive profile basis or its flux-coordinate derivatives.

    $$B_j(u)={n\choose j}u^j(1-u)^{n-j},\qquad u=\psi_N/d$$

    The spline alternative uses six clamped cubic B-splines [1]_. Ordering is
    ascending Bernstein index or spline index; derivatives are in psi_N.

    Parameters
    ----------
    psi_n : array_like
        Finite normalized poloidal flux coordinates [-].
    config : RegularizedProfileConfig
        Basis, domain and knot declaration; scale is unused here [-].
    derivative : int, optional
        Derivative order, zero, one or two [-].

    Returns
    -------
    basis : ndarray
        Coordinates by ordered basis coefficients, with flux derivatives [-].

    Limitations
    -----------
    Extrapolated basis values are mathematical evaluations only. Positive
    partition-of-unity controls guarantee positivity only on the declared domain.

    References
    ----------
    .. [1] C. de Boor, A Practical Guide to Splines, revised edition (2001),
       Chapters IX-X, B-spline representation and differentiation.
    """
    x = np.asarray(psi_n, dtype=float)
    if not np.isfinite(x).all() or derivative not in (0, 1, 2):
        raise ValueError("finite coordinates and derivative 0, 1 or 2 required")
    u = x / config.domain_max
    if config.basis == "bspline6":
        knots = [0.] * 4 + list(config.knot_fractions) + [1.] * 4
        b = BSpline(knots, np.eye(6), 3, extrapolate=True)
        return b(u, nu=derivative) / config.domain_max**derivative
    n = config.coefficient_count - 1
    # Polynomial coefficients give exact derivatives including both endpoints.
    from numpy.polynomial import Polynomial
    cols = []
    for j in range(n + 1):
        p = comb(n, j) * Polynomial([0., 1.])**j * Polynomial([1., -1.])**(n-j)
        cols.append(p.deriv(derivative)(u) / config.domain_max**derivative)
    return np.stack(cols, axis=-1)