image-match/python/py/Lib/site-packages/scipy/_lib/cobyqa/main.py

import warnings

import numpy as np
from scipy.optimize import (
    Bounds,
    LinearConstraint,
    NonlinearConstraint,
    OptimizeResult,
)

from .framework import TrustRegion
from .problem import (
    ObjectiveFunction,
    BoundConstraints,
    LinearConstraints,
    NonlinearConstraints,
    Problem,
)
from .utils import (
    MaxEvalError,
    TargetSuccess,
    CallbackSuccess,
    FeasibleSuccess,
    exact_1d_array,
)
from .settings import (
    ExitStatus,
    Options,
    Constants,
    DEFAULT_OPTIONS,
    DEFAULT_CONSTANTS,
    PRINT_OPTIONS,
)


def minimize(
    fun,
    x0,
    args=(),
    bounds=None,
    constraints=(),
    callback=None,
    options=None,
    **kwargs,
):
    r"""
    Minimize a scalar function using the COBYQA method.

    The Constrained Optimization BY Quadratic Approximations (COBYQA) method is
    a derivative-free optimization method designed to solve general nonlinear
    optimization problems. A complete description of COBYQA is given in [3]_.

    Parameters
    ----------
    fun : {callable, None}
        Objective function to be minimized.

            ``fun(x, *args) -> float``

        where ``x`` is an array with shape (n,) and `args` is a tuple. If `fun`
        is ``None``, the objective function is assumed to be the zero function,
        resulting in a feasibility problem.
    x0 : array_like, shape (n,)
        Initial guess.
    args : tuple, optional
        Extra arguments passed to the objective function.
    bounds : {`scipy.optimize.Bounds`, array_like, shape (n, 2)}, optional
        Bound constraints of the problem. It can be one of the cases below.

        #. An instance of `scipy.optimize.Bounds`. For the time being, the
           argument ``keep_feasible`` is disregarded, and all the constraints
           are considered unrelaxable and will be enforced.
        #. An array with shape (n, 2). The bound constraints for ``x[i]`` are
           ``bounds[i][0] <= x[i] <= bounds[i][1]``. Set ``bounds[i][0]`` to
           :math:`-\infty` if there is no lower bound, and set ``bounds[i][1]``
           to :math:`\infty` if there is no upper bound.

        The COBYQA method always respect the bound constraints.
    constraints : {Constraint, list}, optional
        General constraints of the problem. It can be one of the cases below.

        #. An instance of `scipy.optimize.LinearConstraint`. The argument
           ``keep_feasible`` is disregarded.
        #. An instance of `scipy.optimize.NonlinearConstraint`. The arguments
           ``jac``, ``hess``, ``keep_feasible``, ``finite_diff_rel_step``, and
           ``finite_diff_jac_sparsity`` are disregarded.

        #. A list, each of whose elements are described in the cases above.

    callback : callable, optional
        A callback executed at each objective function evaluation. The method
        terminates if a ``StopIteration`` exception is raised by the callback
        function. Its signature can be one of the following:

            ``callback(intermediate_result)``

        where ``intermediate_result`` is a keyword parameter that contains an
        instance of `scipy.optimize.OptimizeResult`, with attributes ``x``
        and ``fun``, being the point at which the objective function is
        evaluated and the value of the objective function, respectively. The
        name of the parameter must be ``intermediate_result`` for the callback
        to be passed an instance of `scipy.optimize.OptimizeResult`.

        Alternatively, the callback function can have the signature:

            ``callback(xk)``

        where ``xk`` is the point at which the objective function is evaluated.
        Introspection is used to determine which of the signatures to invoke.
    options : dict, optional
        Options passed to the solver. Accepted keys are:

            disp : bool, optional
                Whether to print information about the optimization procedure.
                Default is ``False``.
            maxfev : int, optional
                Maximum number of function evaluations. Default is ``500 * n``.
            maxiter : int, optional
                Maximum number of iterations. Default is ``1000 * n``.
            target : float, optional
                Target on the objective function value. The optimization
                procedure is terminated when the objective function value of a
                feasible point is less than or equal to this target. Default is
                ``-numpy.inf``.
            feasibility_tol : float, optional
                Tolerance on the constraint violation. If the maximum
                constraint violation at a point is less than or equal to this
                tolerance, the point is considered feasible. Default is
                ``numpy.sqrt(numpy.finfo(float).eps)``.
            radius_init : float, optional
                Initial trust-region radius. Typically, this value should be in
                the order of one tenth of the greatest expected change to `x0`.
                Default is ``1.0``.
            radius_final : float, optional
                Final trust-region radius. It should indicate the accuracy
                required in the final values of the variables. Default is
                ``1e-6``.
            nb_points : int, optional
                Number of interpolation points used to build the quadratic
                models of the objective and constraint functions. Default is
                ``2 * n + 1``.
            scale : bool, optional
                Whether to scale the variables according to the bounds. Default
                is ``False``.
            filter_size : int, optional
                Maximum number of points in the filter. The filter is used to
                select the best point returned by the optimization procedure.
                Default is ``sys.maxsize``.
            store_history : bool, optional
                Whether to store the history of the function evaluations.
                Default is ``False``.
            history_size : int, optional
                Maximum number of function evaluations to store in the history.
                Default is ``sys.maxsize``.
            debug : bool, optional
                Whether to perform additional checks during the optimization
                procedure. This option should be used only for debugging
                purposes and is highly discouraged to general users. Default is
                ``False``.

        Other constants (from the keyword arguments) are described below. They
        are not intended to be changed by general users. They should only be
        changed by users with a deep understanding of the algorithm, who want
        to experiment with different settings.

    Returns
    -------
    `scipy.optimize.OptimizeResult`
        Result of the optimization procedure, with the following fields:

            message : str
                Description of the cause of the termination.
            success : bool
                Whether the optimization procedure terminated successfully.
            status : int
                Termination status of the optimization procedure.
            x : `numpy.ndarray`, shape (n,)
                Solution point.
            fun : float
                Objective function value at the solution point.
            maxcv : float
                Maximum constraint violation at the solution point.
            nfev : int
                Number of function evaluations.
            nit : int
                Number of iterations.

        If ``store_history`` is True, the result also has the following fields:

            fun_history : `numpy.ndarray`, shape (nfev,)
                History of the objective function values.
            maxcv_history : `numpy.ndarray`, shape (nfev,)
                History of the maximum constraint violations.

        A description of the termination statuses is given below.

        .. list-table::
            :widths: 25 75
            :header-rows: 1

            * - Exit status
              - Description
            * - 0
              - The lower bound for the trust-region radius has been reached.
            * - 1
              - The target objective function value has been reached.
            * - 2
              - All variables are fixed by the bound constraints.
            * - 3
              - The callback requested to stop the optimization procedure.
            * - 4
              - The feasibility problem received has been solved successfully.
            * - 5
              - The maximum number of function evaluations has been exceeded.
            * - 6
              - The maximum number of iterations has been exceeded.
            * - -1
              - The bound constraints are infeasible.
            * - -2
              - A linear algebra error occurred.

    Other Parameters
    ----------------
    decrease_radius_factor : float, optional
        Factor by which the trust-region radius is reduced when the reduction
        ratio is low or negative. Default is ``0.5``.
    increase_radius_factor : float, optional
        Factor by which the trust-region radius is increased when the reduction
        ratio is large. Default is ``numpy.sqrt(2.0)``.
    increase_radius_threshold : float, optional
        Threshold that controls the increase of the trust-region radius when
        the reduction ratio is large. Default is ``2.0``.
    decrease_radius_threshold : float, optional
        Threshold used to determine whether the trust-region radius should be
        reduced to the resolution. Default is ``1.4``.
    decrease_resolution_factor : float, optional
        Factor by which the resolution is reduced when the current value is far
        from its final value. Default is ``0.1``.
    large_resolution_threshold : float, optional
        Threshold used to determine whether the resolution is far from its
        final value. Default is ``250.0``.
    moderate_resolution_threshold : float, optional
        Threshold used to determine whether the resolution is close to its
        final value. Default is ``16.0``.
    low_ratio : float, optional
        Threshold used to determine whether the reduction ratio is low. Default
        is ``0.1``.
    high_ratio : float, optional
        Threshold used to determine whether the reduction ratio is high.
        Default is ``0.7``.
    very_low_ratio : float, optional
        Threshold used to determine whether the reduction ratio is very low.
        This is used to determine whether the models should be reset. Default
        is ``0.01``.
    penalty_increase_threshold : float, optional
        Threshold used to determine whether the penalty parameter should be
        increased. Default is ``1.5``.
    penalty_increase_factor : float, optional
        Factor by which the penalty parameter is increased. Default is ``2.0``.
    short_step_threshold : float, optional
        Factor used to determine whether the trial step is too short. Default
        is ``0.5``.
    low_radius_factor : float, optional
        Factor used to determine which interpolation point should be removed
        from the interpolation set at each iteration. Default is ``0.1``.
    byrd_omojokun_factor : float, optional
        Factor by which the trust-region radius is reduced for the computations
        of the normal step in the Byrd-Omojokun composite-step approach.
        Default is ``0.8``.
    threshold_ratio_constraints : float, optional
        Threshold used to determine which constraints should be taken into
        account when decreasing the penalty parameter. Default is ``2.0``.
    large_shift_factor : float, optional
        Factor used to determine whether the point around which the quadratic
        models are built should be updated. Default is ``10.0``.
    large_gradient_factor : float, optional
        Factor used to determine whether the models should be reset. Default is
        ``10.0``.
    resolution_factor : float, optional
        Factor by which the resolution is decreased. Default is ``2.0``.
    improve_tcg : bool, optional
        Whether to improve the steps computed by the truncated conjugate
        gradient method when the trust-region boundary is reached. Default is
        ``True``.

    References
    ----------
    .. [1] J. Nocedal and S. J. Wright. *Numerical Optimization*. Springer Ser.
       Oper. Res. Financ. Eng. Springer, New York, NY, USA, second edition,
       2006. `doi:10.1007/978-0-387-40065-5
       <https://doi.org/10.1007/978-0-387-40065-5>`_.
    .. [2] M. J. D. Powell. A direct search optimization method that models the
       objective and constraint functions by linear interpolation. In S. Gomez
       and J.-P. Hennart, editors, *Advances in Optimization and Numerical
       Analysis*, volume 275 of Math. Appl., pages 51--67. Springer, Dordrecht,
       Netherlands, 1994. `doi:10.1007/978-94-015-8330-5_4
       <https://doi.org/10.1007/978-94-015-8330-5_4>`_.
    .. [3] T. M. Ragonneau. *Model-Based Derivative-Free Optimization Methods
       and Software*. PhD thesis, Department of Applied Mathematics, The Hong
       Kong Polytechnic University, Hong Kong, China, 2022. URL:
       https://theses.lib.polyu.edu.hk/handle/200/12294.

    Examples
    --------
    To demonstrate how to use `minimize`, we first minimize the Rosenbrock
    function implemented in `scipy.optimize` in an unconstrained setting.

    .. testsetup::

        import numpy as np
        np.set_printoptions(precision=3, suppress=True)

    >>> from cobyqa import minimize
    >>> from scipy.optimize import rosen

    To solve the problem using COBYQA, run:

    >>> x0 = [1.3, 0.7, 0.8, 1.9, 1.2]
    >>> res = minimize(rosen, x0)
    >>> res.x
    array([1., 1., 1., 1., 1.])

    To see how bound and constraints are handled using `minimize`, we solve
    Example 16.4 of [1]_, defined as

    .. math::

        \begin{aligned}
            \min_{x \in \mathbb{R}^2}   & \quad (x_1 - 1)^2 + (x_2 - 2.5)^2\\
            \text{s.t.}                 & \quad -x_1 + 2x_2 \le 2,\\
                                        & \quad x_1 + 2x_2 \le 6,\\
                                        & \quad x_1 - 2x_2 \le 2,\\
                                        & \quad x_1 \ge 0,\\
                                        & \quad x_2 \ge 0.
        \end{aligned}

    >>> import numpy as np
    >>> from scipy.optimize import Bounds, LinearConstraint

    Its objective function can be implemented as:

    >>> def fun(x):
    ...     return (x[0] - 1.0)**2 + (x[1] - 2.5)**2

    This problem can be solved using `minimize` as:

    >>> x0 = [2.0, 0.0]
    >>> bounds = Bounds([0.0, 0.0], np.inf)
    >>> constraints = LinearConstraint([
    ...     [-1.0, 2.0],
    ...     [1.0, 2.0],
    ...     [1.0, -2.0],
    ... ], -np.inf, [2.0, 6.0, 2.0])
    >>> res = minimize(fun, x0, bounds=bounds, constraints=constraints)
    >>> res.x
    array([1.4, 1.7])

    To see how nonlinear constraints are handled, we solve Problem (F) of [2]_,
    defined as

    .. math::

        \begin{aligned}
            \min_{x \in \mathbb{R}^2}   & \quad -x_1 - x_2\\
            \text{s.t.}                 & \quad x_1^2 - x_2 \le 0,\\
                                        & \quad x_1^2 + x_2^2 \le 1.
        \end{aligned}

    >>> from scipy.optimize import NonlinearConstraint

    Its objective and constraint functions can be implemented as:

    >>> def fun(x):
    ...     return -x[0] - x[1]
    >>>
    >>> def cub(x):
    ...     return [x[0]**2 - x[1], x[0]**2 + x[1]**2]

    This problem can be solved using `minimize` as:

    >>> x0 = [1.0, 1.0]
    >>> constraints = NonlinearConstraint(cub, -np.inf, [0.0, 1.0])
    >>> res = minimize(fun, x0, constraints=constraints)
    >>> res.x
    array([0.707, 0.707])

    Finally, to see how to supply linear and nonlinear constraints
    simultaneously, we solve Problem (G) of [2]_, defined as

    .. math::

        \begin{aligned}
            \min_{x \in \mathbb{R}^3}   & \quad x_3\\
            \text{s.t.}                 & \quad 5x_1 - x_2 + x_3 \ge 0,\\
                                        & \quad -5x_1 - x_2 + x_3 \ge 0,\\
                                        & \quad x_1^2 + x_2^2 + 4x_2 \le x_3.
        \end{aligned}

    Its objective and nonlinear constraint functions can be implemented as:

    >>> def fun(x):
    ...     return x[2]
    >>>
    >>> def cub(x):
    ...     return x[0]**2 + x[1]**2 + 4.0*x[1] - x[2]

    This problem can be solved using `minimize` as:

    >>> x0 = [1.0, 1.0, 1.0]
    >>> constraints = [
    ...     LinearConstraint(
    ...         [[5.0, -1.0, 1.0], [-5.0, -1.0, 1.0]],
    ...         [0.0, 0.0],
    ...         np.inf,
    ...     ),
    ...     NonlinearConstraint(cub, -np.inf, 0.0),
    ... ]
    >>> res = minimize(fun, x0, constraints=constraints)
    >>> res.x
    array([ 0., -3., -3.])
    """
    # Get basic options that are needed for the initialization.
    if options is None:
        options = {}
    else:
        options = dict(options)
    verbose = options.get(Options.VERBOSE, DEFAULT_OPTIONS[Options.VERBOSE])
    verbose = bool(verbose)
    feasibility_tol = options.get(
        Options.FEASIBILITY_TOL,
        DEFAULT_OPTIONS[Options.FEASIBILITY_TOL],
    )
    feasibility_tol = float(feasibility_tol)
    scale = options.get(Options.SCALE, DEFAULT_OPTIONS[Options.SCALE])
    scale = bool(scale)
    store_history = options.get(
        Options.STORE_HISTORY,
        DEFAULT_OPTIONS[Options.STORE_HISTORY],
    )
    store_history = bool(store_history)
    if Options.HISTORY_SIZE in options and options[Options.HISTORY_SIZE] <= 0:
        raise ValueError("The size of the history must be positive.")
    history_size = options.get(
        Options.HISTORY_SIZE,
        DEFAULT_OPTIONS[Options.HISTORY_SIZE],
    )
    history_size = int(history_size)
    if Options.FILTER_SIZE in options and options[Options.FILTER_SIZE] <= 0:
        raise ValueError("The size of the filter must be positive.")
    filter_size = options.get(
        Options.FILTER_SIZE,
        DEFAULT_OPTIONS[Options.FILTER_SIZE],
    )
    filter_size = int(filter_size)
    debug = options.get(Options.DEBUG, DEFAULT_OPTIONS[Options.DEBUG])
    debug = bool(debug)

    # Initialize the objective function.
    if not isinstance(args, tuple):
        args = (args,)
    obj = ObjectiveFunction(fun, verbose, debug, *args)

    # Initialize the bound constraints.
    if not hasattr(x0, "__len__"):
        x0 = [x0]
    n_orig = len(x0)
    bounds = BoundConstraints(_get_bounds(bounds, n_orig))

    # Initialize the constraints.
    linear_constraints, nonlinear_constraints = _get_constraints(constraints)
    linear = LinearConstraints(linear_constraints, n_orig, debug)
    nonlinear = NonlinearConstraints(nonlinear_constraints, verbose, debug)

    # Initialize the problem (and remove the fixed variables).
    pb = Problem(
        obj,
        x0,
        bounds,
        linear,
        nonlinear,
        callback,
        feasibility_tol,
        scale,
        store_history,
        history_size,
        filter_size,
        debug,
    )

    # Set the default options.
    _set_default_options(options, pb.n)
    constants = _set_default_constants(**kwargs)

    # Initialize the models and skip the computations whenever possible.
    if not pb.bounds.is_feasible:
        # The bound constraints are infeasible.
        return _build_result(
            pb,
            0.0,
            False,
            ExitStatus.INFEASIBLE_ERROR,
            0,
            options,
        )
    elif pb.n == 0:
        # All variables are fixed by the bound constraints.
        return _build_result(
            pb,
            0.0,
            True,
            ExitStatus.FIXED_SUCCESS,
            0,
            options,
        )
    if verbose:
        print("Starting the optimization procedure.")
        print(f"Initial trust-region radius: {options[Options.RHOBEG]}.")
        print(f"Final trust-region radius: {options[Options.RHOEND]}.")
        print(
            f"Maximum number of function evaluations: "
            f"{options[Options.MAX_EVAL]}."
        )
        print(f"Maximum number of iterations: {options[Options.MAX_ITER]}.")
        print()
    try:
        framework = TrustRegion(pb, options, constants)
    except TargetSuccess:
        # The target on the objective function value has been reached
        return _build_result(
            pb,
            0.0,
            True,
            ExitStatus.TARGET_SUCCESS,
            0,
            options,
        )
    except CallbackSuccess:
        # The callback raised a StopIteration exception.
        return _build_result(
            pb,
            0.0,
            True,
            ExitStatus.CALLBACK_SUCCESS,
            0,
            options,
        )
    except FeasibleSuccess:
        # The feasibility problem has been solved successfully.
        return _build_result(
            pb,
            0.0,
            True,
            ExitStatus.FEASIBLE_SUCCESS,
            0,
            options,
        )
    except MaxEvalError:
        # The maximum number of function evaluations has been exceeded.
        return _build_result(
            pb,
            0.0,
            False,
            ExitStatus.MAX_ITER_WARNING,
            0,
            options,
        )
    except np.linalg.LinAlgError:
        # The construction of the initial interpolation set failed.
        return _build_result(
            pb,
            0.0,
            False,
            ExitStatus.LINALG_ERROR,
            0,
            options,
        )

    # Start the optimization procedure.
    success = False
    n_iter = 0
    k_new = None
    n_short_steps = 0
    n_very_short_steps = 0
    n_alt_models = 0
    while True:
        # Stop the optimization procedure if the maximum number of iterations
        # has been exceeded. We do not write the main loop as a for loop
        # because we want to access the number of iterations outside the loop.
        if n_iter >= options[Options.MAX_ITER]:
            status = ExitStatus.MAX_ITER_WARNING
            break
        n_iter += 1

        # Update the point around which the quadratic models are built.
        if (
            np.linalg.norm(
                framework.x_best - framework.models.interpolation.x_base
            )
            >= constants[Constants.LARGE_SHIFT_FACTOR] * framework.radius
        ):
            framework.shift_x_base(options)

        # Evaluate the trial step.
        radius_save = framework.radius
        normal_step, tangential_step = framework.get_trust_region_step(options)
        step = normal_step + tangential_step
        s_norm = np.linalg.norm(step)

        # If the trial step is too short, we do not attempt to evaluate the
        # objective and constraint functions. Instead, we reduce the
        # trust-region radius and check whether the resolution should be
        # enhanced and whether the geometry of the interpolation set should be
        # improved. Otherwise, we entertain a classical iteration. The
        # criterion for performing an exceptional jump is taken from NEWUOA.
        if (
            s_norm
            <= constants[Constants.SHORT_STEP_THRESHOLD] * framework.resolution
        ):
            framework.radius *= constants[Constants.DECREASE_RESOLUTION_FACTOR]
            if radius_save > framework.resolution:
                n_short_steps = 0
                n_very_short_steps = 0
            else:
                n_short_steps += 1
                n_very_short_steps += 1
                if s_norm > 0.1 * framework.resolution:
                    n_very_short_steps = 0
            enhance_resolution = n_short_steps >= 5 or n_very_short_steps >= 3
            if enhance_resolution:
                n_short_steps = 0
                n_very_short_steps = 0
                improve_geometry = False
            else:
                try:
                    k_new, dist_new = framework.get_index_to_remove()
                except np.linalg.LinAlgError:
                    status = ExitStatus.LINALG_ERROR
                    break
                improve_geometry = dist_new > max(
                    framework.radius,
                    constants[Constants.RESOLUTION_FACTOR]
                    * framework.resolution,
                )
        else:
            # Increase the penalty parameter if necessary.
            same_best_point = framework.increase_penalty(step)
            if same_best_point:
                # Evaluate the objective and constraint functions.
                try:
                    fun_val, cub_val, ceq_val = _eval(
                        pb,
                        framework,
                        step,
                        options,
                    )
                except TargetSuccess:
                    status = ExitStatus.TARGET_SUCCESS
                    success = True
                    break
                except FeasibleSuccess:
                    status = ExitStatus.FEASIBLE_SUCCESS
                    success = True
                    break
                except CallbackSuccess:
                    status = ExitStatus.CALLBACK_SUCCESS
                    success = True
                    break
                except MaxEvalError:
                    status = ExitStatus.MAX_EVAL_WARNING
                    break

                # Perform a second-order correction step if necessary.
                merit_old = framework.merit(
                    framework.x_best,
                    framework.fun_best,
                    framework.cub_best,
                    framework.ceq_best,
                )
                merit_new = framework.merit(
                    framework.x_best + step, fun_val, cub_val, ceq_val
                )
                if (
                    pb.type == "nonlinearly constrained"
                    and merit_new > merit_old
                    and np.linalg.norm(normal_step)
                    > constants[Constants.BYRD_OMOJOKUN_FACTOR] ** 2.0
                    * framework.radius
                ):
                    soc_step = framework.get_second_order_correction_step(
                        step, options
                    )
                    if np.linalg.norm(soc_step) > 0.0:
                        step += soc_step

                        # Evaluate the objective and constraint functions.
                        try:
                            fun_val, cub_val, ceq_val = _eval(
                                pb,
                                framework,
                                step,
                                options,
                            )
                        except TargetSuccess:
                            status = ExitStatus.TARGET_SUCCESS
                            success = True
                            break
                        except FeasibleSuccess:
                            status = ExitStatus.FEASIBLE_SUCCESS
                            success = True
                            break
                        except CallbackSuccess:
                            status = ExitStatus.CALLBACK_SUCCESS
                            success = True
                            break
                        except MaxEvalError:
                            status = ExitStatus.MAX_EVAL_WARNING
                            break

                # Calculate the reduction ratio.
                ratio = framework.get_reduction_ratio(
                    step,
                    fun_val,
                    cub_val,
                    ceq_val,
                )

                # Choose an interpolation point to remove.
                try:
                    k_new = framework.get_index_to_remove(
                        framework.x_best + step
                    )[0]
                except np.linalg.LinAlgError:
                    status = ExitStatus.LINALG_ERROR
                    break

                # Update the interpolation set.
                try:
                    ill_conditioned = framework.models.update_interpolation(
                        k_new, framework.x_best + step, fun_val, cub_val,
                        ceq_val
                    )
                except np.linalg.LinAlgError:
                    status = ExitStatus.LINALG_ERROR
                    break
                framework.set_best_index()

                # Update the trust-region radius.
                framework.update_radius(step, ratio)

                # Attempt to replace the models by the alternative ones.
                if framework.radius <= framework.resolution:
                    if ratio >= constants[Constants.VERY_LOW_RATIO]:
                        n_alt_models = 0
                    else:
                        n_alt_models += 1
                        grad = framework.models.fun_grad(framework.x_best)
                        try:
                            grad_alt = framework.models.fun_alt_grad(
                                framework.x_best
                            )
                        except np.linalg.LinAlgError:
                            status = ExitStatus.LINALG_ERROR
                            break
                        if np.linalg.norm(grad) < constants[
                            Constants.LARGE_GRADIENT_FACTOR
                        ] * np.linalg.norm(grad_alt):
                            n_alt_models = 0
                        if n_alt_models >= 3:
                            try:
                                framework.models.reset_models()
                            except np.linalg.LinAlgError:
                                status = ExitStatus.LINALG_ERROR
                                break
                            n_alt_models = 0

                # Update the Lagrange multipliers.
                framework.set_multipliers(framework.x_best + step)

                # Check whether the resolution should be enhanced.
                try:
                    k_new, dist_new = framework.get_index_to_remove()
                except np.linalg.LinAlgError:
                    status = ExitStatus.LINALG_ERROR
                    break
                improve_geometry = (
                    ill_conditioned
                    or ratio <= constants[Constants.LOW_RATIO]
                    and dist_new
                    > max(
                        framework.radius,
                        constants[Constants.RESOLUTION_FACTOR]
                        * framework.resolution,
                    )
                )
                enhance_resolution = (
                    radius_save <= framework.resolution
                    and ratio <= constants[Constants.LOW_RATIO]
                    and not improve_geometry
                )
            else:
                # When increasing the penalty parameter, the best point so far
                # may change. In this case, we restart the iteration.
                enhance_resolution = False
                improve_geometry = False

        # Reduce the resolution if necessary.
        if enhance_resolution:
            if framework.resolution <= options[Options.RHOEND]:
                success = True
                status = ExitStatus.RADIUS_SUCCESS
                break
            framework.enhance_resolution(options)
            framework.decrease_penalty()

            if verbose:
                maxcv_val = pb.maxcv(
                    framework.x_best, framework.cub_best, framework.ceq_best
                )
                _print_step(
                    f"New trust-region radius: {framework.resolution}",
                    pb,
                    pb.build_x(framework.x_best),
                    framework.fun_best,
                    maxcv_val,
                    pb.n_eval,
                    n_iter,
                )
                print()

        # Improve the geometry of the interpolation set if necessary.
        if improve_geometry:
            try:
                step = framework.get_geometry_step(k_new, options)
            except np.linalg.LinAlgError:
                status = ExitStatus.LINALG_ERROR
                break

            # Evaluate the objective and constraint functions.
            try:
                fun_val, cub_val, ceq_val = _eval(pb, framework, step, options)
            except TargetSuccess:
                status = ExitStatus.TARGET_SUCCESS
                success = True
                break
            except FeasibleSuccess:
                status = ExitStatus.FEASIBLE_SUCCESS
                success = True
                break
            except CallbackSuccess:
                status = ExitStatus.CALLBACK_SUCCESS
                success = True
                break
            except MaxEvalError:
                status = ExitStatus.MAX_EVAL_WARNING
                break

            # Update the interpolation set.
            try:
                framework.models.update_interpolation(
                    k_new,
                    framework.x_best + step,
                    fun_val,
                    cub_val,
                    ceq_val,
                )
            except np.linalg.LinAlgError:
                status = ExitStatus.LINALG_ERROR
                break
            framework.set_best_index()

    return _build_result(
        pb,
        framework.penalty,
        success,
        status,
        n_iter,
        options,
    )


def _get_bounds(bounds, n):
    """
    Uniformize the bounds.
    """
    if bounds is None:
        return Bounds(np.full(n, -np.inf), np.full(n, np.inf))
    elif isinstance(bounds, Bounds):
        if bounds.lb.shape != (n,) or bounds.ub.shape != (n,):
            raise ValueError(f"The bounds must have {n} elements.")
        return Bounds(bounds.lb, bounds.ub)
    elif hasattr(bounds, "__len__"):
        bounds = np.asarray(bounds)
        if bounds.shape != (n, 2):
            raise ValueError(
                "The shape of the bounds is not compatible with "
                "the number of variables."
            )
        return Bounds(bounds[:, 0], bounds[:, 1])
    else:
        raise TypeError(
            "The bounds must be an instance of "
            "scipy.optimize.Bounds or an array-like object."
        )


def _get_constraints(constraints):
    """
    Extract the linear and nonlinear constraints.
    """
    if isinstance(constraints, dict) or not hasattr(constraints, "__len__"):
        constraints = (constraints,)

    # Extract the linear and nonlinear constraints.
    linear_constraints = []
    nonlinear_constraints = []
    for constraint in constraints:
        if isinstance(constraint, LinearConstraint):
            lb = exact_1d_array(
                constraint.lb,
                "The lower bound of the linear constraints must be a vector.",
            )
            ub = exact_1d_array(
                constraint.ub,
                "The upper bound of the linear constraints must be a vector.",
            )
            linear_constraints.append(
                LinearConstraint(
                    constraint.A,
                    *np.broadcast_arrays(lb, ub),
                )
            )
        elif isinstance(constraint, NonlinearConstraint):
            lb = exact_1d_array(
                constraint.lb,
                "The lower bound of the "
                "nonlinear constraints must be a "
                "vector.",
            )
            ub = exact_1d_array(
                constraint.ub,
                "The upper bound of the "
                "nonlinear constraints must be a "
                "vector.",
            )
            nonlinear_constraints.append(
                NonlinearConstraint(
                    constraint.fun,
                    *np.broadcast_arrays(lb, ub),
                )
            )
        elif isinstance(constraint, dict):
            if "type" not in constraint or constraint["type"] not in (
                "eq",
                "ineq",
            ):
                raise ValueError('The constraint type must be "eq" or "ineq".')
            if "fun" not in constraint or not callable(constraint["fun"]):
                raise ValueError("The constraint function must be callable.")
            nonlinear_constraints.append(
                {
                    "fun": constraint["fun"],
                    "type": constraint["type"],
                    "args": constraint.get("args", ()),
                }
            )
        else:
            raise TypeError(
                "The constraints must be instances of "
                "scipy.optimize.LinearConstraint, "
                "scipy.optimize.NonlinearConstraint, or dict."
            )
    return linear_constraints, nonlinear_constraints


def _set_default_options(options, n):
    """
    Set the default options.
    """
    if Options.RHOBEG in options and options[Options.RHOBEG] <= 0.0:
        raise ValueError("The initial trust-region radius must be positive.")
    if Options.RHOEND in options and options[Options.RHOEND] < 0.0:
        raise ValueError("The final trust-region radius must be nonnegative.")
    if Options.RHOBEG in options and Options.RHOEND in options:
        if options[Options.RHOBEG] < options[Options.RHOEND]:
            raise ValueError(
                "The initial trust-region radius must be greater "
                "than or equal to the final trust-region radius."
            )
    elif Options.RHOBEG in options:
        options[Options.RHOEND.value] = np.min(
            [
                DEFAULT_OPTIONS[Options.RHOEND],
                options[Options.RHOBEG],
            ]
        )
    elif Options.RHOEND in options:
        options[Options.RHOBEG.value] = np.max(
            [
                DEFAULT_OPTIONS[Options.RHOBEG],
                options[Options.RHOEND],
            ]
        )
    else:
        options[Options.RHOBEG.value] = DEFAULT_OPTIONS[Options.RHOBEG]
        options[Options.RHOEND.value] = DEFAULT_OPTIONS[Options.RHOEND]
    options[Options.RHOBEG.value] = float(options[Options.RHOBEG])
    options[Options.RHOEND.value] = float(options[Options.RHOEND])
    if Options.NPT in options and options[Options.NPT] <= 0:
        raise ValueError("The number of interpolation points must be "
                         "positive.")
    if (
        Options.NPT in options
        and options[Options.NPT] > ((n + 1) * (n + 2)) // 2
    ):
        raise ValueError(
            f"The number of interpolation points must be at most "
            f"{((n + 1) * (n + 2)) // 2}."
        )
    options.setdefault(Options.NPT.value, DEFAULT_OPTIONS[Options.NPT](n))
    options[Options.NPT.value] = int(options[Options.NPT])
    if Options.MAX_EVAL in options and options[Options.MAX_EVAL] <= 0:
        raise ValueError(
            "The maximum number of function evaluations must be positive."
        )
    options.setdefault(
        Options.MAX_EVAL.value,
        np.max(
            [
                DEFAULT_OPTIONS[Options.MAX_EVAL](n),
                options[Options.NPT] + 1,
            ]
        ),
    )
    options[Options.MAX_EVAL.value] = int(options[Options.MAX_EVAL])
    if Options.MAX_ITER in options and options[Options.MAX_ITER] <= 0:
        raise ValueError("The maximum number of iterations must be positive.")
    options.setdefault(
        Options.MAX_ITER.value,
        DEFAULT_OPTIONS[Options.MAX_ITER](n),
    )
    options[Options.MAX_ITER.value] = int(options[Options.MAX_ITER])
    options.setdefault(Options.TARGET.value, DEFAULT_OPTIONS[Options.TARGET])
    options[Options.TARGET.value] = float(options[Options.TARGET])
    options.setdefault(
        Options.FEASIBILITY_TOL.value,
        DEFAULT_OPTIONS[Options.FEASIBILITY_TOL],
    )
    options[Options.FEASIBILITY_TOL.value] = float(
        options[Options.FEASIBILITY_TOL]
    )
    options.setdefault(Options.VERBOSE.value, DEFAULT_OPTIONS[Options.VERBOSE])
    options[Options.VERBOSE.value] = bool(options[Options.VERBOSE])
    options.setdefault(Options.SCALE.value, DEFAULT_OPTIONS[Options.SCALE])
    options[Options.SCALE.value] = bool(options[Options.SCALE])
    options.setdefault(
        Options.FILTER_SIZE.value,
        DEFAULT_OPTIONS[Options.FILTER_SIZE],
    )
    options[Options.FILTER_SIZE.value] = int(options[Options.FILTER_SIZE])
    options.setdefault(
        Options.STORE_HISTORY.value,
        DEFAULT_OPTIONS[Options.STORE_HISTORY],
    )
    options[Options.STORE_HISTORY.value] = bool(options[Options.STORE_HISTORY])
    options.setdefault(
        Options.HISTORY_SIZE.value,
        DEFAULT_OPTIONS[Options.HISTORY_SIZE],
    )
    options[Options.HISTORY_SIZE.value] = int(options[Options.HISTORY_SIZE])
    options.setdefault(Options.DEBUG.value, DEFAULT_OPTIONS[Options.DEBUG])
    options[Options.DEBUG.value] = bool(options[Options.DEBUG])

    # Check whether they are any unknown options.
    for key in options:
        if key not in Options.__members__.values():
            warnings.warn(f"Unknown option: {key}.", RuntimeWarning, 3)


def _set_default_constants(**kwargs):
    """
    Set the default constants.
    """
    constants = dict(kwargs)
    constants.setdefault(
        Constants.DECREASE_RADIUS_FACTOR.value,
        DEFAULT_CONSTANTS[Constants.DECREASE_RADIUS_FACTOR],
    )
    constants[Constants.DECREASE_RADIUS_FACTOR.value] = float(
        constants[Constants.DECREASE_RADIUS_FACTOR]
    )
    if (
        constants[Constants.DECREASE_RADIUS_FACTOR] <= 0.0
        or constants[Constants.DECREASE_RADIUS_FACTOR] >= 1.0
    ):
        raise ValueError(
            "The constant decrease_radius_factor must be in the interval "
            "(0, 1)."
        )
    constants.setdefault(
        Constants.INCREASE_RADIUS_THRESHOLD.value,
        DEFAULT_CONSTANTS[Constants.INCREASE_RADIUS_THRESHOLD],
    )
    constants[Constants.INCREASE_RADIUS_THRESHOLD.value] = float(
        constants[Constants.INCREASE_RADIUS_THRESHOLD]
    )
    if constants[Constants.INCREASE_RADIUS_THRESHOLD] <= 1.0:
        raise ValueError(
            "The constant increase_radius_threshold must be greater than 1."
        )
    if (
        Constants.INCREASE_RADIUS_FACTOR in constants
        and constants[Constants.INCREASE_RADIUS_FACTOR] <= 1.0
    ):
        raise ValueError(
            "The constant increase_radius_factor must be greater than 1."
        )
    if (
        Constants.DECREASE_RADIUS_THRESHOLD in constants
        and constants[Constants.DECREASE_RADIUS_THRESHOLD] <= 1.0
    ):
        raise ValueError(
            "The constant decrease_radius_threshold must be greater than 1."
        )
    if (
        Constants.INCREASE_RADIUS_FACTOR in constants
        and Constants.DECREASE_RADIUS_THRESHOLD in constants
    ):
        if (
            constants[Constants.DECREASE_RADIUS_THRESHOLD]
            >= constants[Constants.INCREASE_RADIUS_FACTOR]
        ):
            raise ValueError(
                "The constant decrease_radius_threshold must be "
                "less than increase_radius_factor."
            )
    elif Constants.INCREASE_RADIUS_FACTOR in constants:
        constants[Constants.DECREASE_RADIUS_THRESHOLD.value] = np.min(
            [
                DEFAULT_CONSTANTS[Constants.DECREASE_RADIUS_THRESHOLD],
                0.5 * (1.0 + constants[Constants.INCREASE_RADIUS_FACTOR]),
            ]
        )
    elif Constants.DECREASE_RADIUS_THRESHOLD in constants:
        constants[Constants.INCREASE_RADIUS_FACTOR.value] = np.max(
            [
                DEFAULT_CONSTANTS[Constants.INCREASE_RADIUS_FACTOR],
                2.0 * constants[Constants.DECREASE_RADIUS_THRESHOLD],
            ]
        )
    else:
        constants[Constants.INCREASE_RADIUS_FACTOR.value] = DEFAULT_CONSTANTS[
            Constants.INCREASE_RADIUS_FACTOR
        ]
        constants[Constants.DECREASE_RADIUS_THRESHOLD.value] = (
            DEFAULT_CONSTANTS[Constants.DECREASE_RADIUS_THRESHOLD])
    constants.setdefault(
        Constants.DECREASE_RESOLUTION_FACTOR.value,
        DEFAULT_CONSTANTS[Constants.DECREASE_RESOLUTION_FACTOR],
    )
    constants[Constants.DECREASE_RESOLUTION_FACTOR.value] = float(
        constants[Constants.DECREASE_RESOLUTION_FACTOR]
    )
    if (
        constants[Constants.DECREASE_RESOLUTION_FACTOR] <= 0.0
        or constants[Constants.DECREASE_RESOLUTION_FACTOR] >= 1.0
    ):
        raise ValueError(
            "The constant decrease_resolution_factor must be in the interval "
            "(0, 1)."
        )
    if (
        Constants.LARGE_RESOLUTION_THRESHOLD in constants
        and constants[Constants.LARGE_RESOLUTION_THRESHOLD] <= 1.0
    ):
        raise ValueError(
            "The constant large_resolution_threshold must be greater than 1."
        )
    if (
        Constants.MODERATE_RESOLUTION_THRESHOLD in constants
        and constants[Constants.MODERATE_RESOLUTION_THRESHOLD] <= 1.0
    ):
        raise ValueError(
            "The constant moderate_resolution_threshold must be greater than "
            "1."
        )
    if (
        Constants.LARGE_RESOLUTION_THRESHOLD in constants
        and Constants.MODERATE_RESOLUTION_THRESHOLD in constants
    ):
        if (
            constants[Constants.MODERATE_RESOLUTION_THRESHOLD]
            > constants[Constants.LARGE_RESOLUTION_THRESHOLD]
        ):
            raise ValueError(
                "The constant moderate_resolution_threshold "
                "must be at most large_resolution_threshold."
            )
    elif Constants.LARGE_RESOLUTION_THRESHOLD in constants:
        constants[Constants.MODERATE_RESOLUTION_THRESHOLD.value] = np.min(
            [
                DEFAULT_CONSTANTS[Constants.MODERATE_RESOLUTION_THRESHOLD],
                constants[Constants.LARGE_RESOLUTION_THRESHOLD],
            ]
        )
    elif Constants.MODERATE_RESOLUTION_THRESHOLD in constants:
        constants[Constants.LARGE_RESOLUTION_THRESHOLD.value] = np.max(
            [
                DEFAULT_CONSTANTS[Constants.LARGE_RESOLUTION_THRESHOLD],
                constants[Constants.MODERATE_RESOLUTION_THRESHOLD],
            ]
        )
    else:
        constants[Constants.LARGE_RESOLUTION_THRESHOLD.value] = (
            DEFAULT_CONSTANTS[Constants.LARGE_RESOLUTION_THRESHOLD]
        )
        constants[Constants.MODERATE_RESOLUTION_THRESHOLD.value] = (
            DEFAULT_CONSTANTS[Constants.MODERATE_RESOLUTION_THRESHOLD]
        )
    if Constants.LOW_RATIO in constants and (
        constants[Constants.LOW_RATIO] <= 0.0
        or constants[Constants.LOW_RATIO] >= 1.0
    ):
        raise ValueError(
            "The constant low_ratio must be in the interval (0, 1)."
        )
    if Constants.HIGH_RATIO in constants and (
        constants[Constants.HIGH_RATIO] <= 0.0
        or constants[Constants.HIGH_RATIO] >= 1.0
    ):
        raise ValueError(
            "The constant high_ratio must be in the interval (0, 1)."
        )
    if Constants.LOW_RATIO in constants and Constants.HIGH_RATIO in constants:
        if constants[Constants.LOW_RATIO] > constants[Constants.HIGH_RATIO]:
            raise ValueError(
                "The constant low_ratio must be at most high_ratio."
            )
    elif Constants.LOW_RATIO in constants:
        constants[Constants.HIGH_RATIO.value] = np.max(
            [
                DEFAULT_CONSTANTS[Constants.HIGH_RATIO],
                constants[Constants.LOW_RATIO],
            ]
        )
    elif Constants.HIGH_RATIO in constants:
        constants[Constants.LOW_RATIO.value] = np.min(
            [
                DEFAULT_CONSTANTS[Constants.LOW_RATIO],
                constants[Constants.HIGH_RATIO],
            ]
        )
    else:
        constants[Constants.LOW_RATIO.value] = DEFAULT_CONSTANTS[
            Constants.LOW_RATIO
        ]
        constants[Constants.HIGH_RATIO.value] = DEFAULT_CONSTANTS[
            Constants.HIGH_RATIO
        ]
    constants.setdefault(
        Constants.VERY_LOW_RATIO.value,
        DEFAULT_CONSTANTS[Constants.VERY_LOW_RATIO],
    )
    constants[Constants.VERY_LOW_RATIO.value] = float(
        constants[Constants.VERY_LOW_RATIO]
    )
    if (
        constants[Constants.VERY_LOW_RATIO] <= 0.0
        or constants[Constants.VERY_LOW_RATIO] >= 1.0
    ):
        raise ValueError(
            "The constant very_low_ratio must be in the interval (0, 1)."
        )
    if (
        Constants.PENALTY_INCREASE_THRESHOLD in constants
        and constants[Constants.PENALTY_INCREASE_THRESHOLD] < 1.0
    ):
        raise ValueError(
            "The constant penalty_increase_threshold must be "
            "greater than or equal to 1."
        )
    if (
        Constants.PENALTY_INCREASE_FACTOR in constants
        and constants[Constants.PENALTY_INCREASE_FACTOR] <= 1.0
    ):
        raise ValueError(
            "The constant penalty_increase_factor must be greater than 1."
        )
    if (
        Constants.PENALTY_INCREASE_THRESHOLD in constants
        and Constants.PENALTY_INCREASE_FACTOR in constants
    ):
        if (
            constants[Constants.PENALTY_INCREASE_FACTOR]
            < constants[Constants.PENALTY_INCREASE_THRESHOLD]
        ):
            raise ValueError(
                "The constant penalty_increase_factor must be "
                "greater than or equal to "
                "penalty_increase_threshold."
            )
    elif Constants.PENALTY_INCREASE_THRESHOLD in constants:
        constants[Constants.PENALTY_INCREASE_FACTOR.value] = np.max(
            [
                DEFAULT_CONSTANTS[Constants.PENALTY_INCREASE_FACTOR],
                constants[Constants.PENALTY_INCREASE_THRESHOLD],
            ]
        )
    elif Constants.PENALTY_INCREASE_FACTOR in constants:
        constants[Constants.PENALTY_INCREASE_THRESHOLD.value] = np.min(
            [
                DEFAULT_CONSTANTS[Constants.PENALTY_INCREASE_THRESHOLD],
                constants[Constants.PENALTY_INCREASE_FACTOR],
            ]
        )
    else:
        constants[Constants.PENALTY_INCREASE_THRESHOLD.value] = (
            DEFAULT_CONSTANTS[Constants.PENALTY_INCREASE_THRESHOLD]
        )
        constants[Constants.PENALTY_INCREASE_FACTOR.value] = DEFAULT_CONSTANTS[
            Constants.PENALTY_INCREASE_FACTOR
        ]
    constants.setdefault(
        Constants.SHORT_STEP_THRESHOLD.value,
        DEFAULT_CONSTANTS[Constants.SHORT_STEP_THRESHOLD],
    )
    constants[Constants.SHORT_STEP_THRESHOLD.value] = float(
        constants[Constants.SHORT_STEP_THRESHOLD]
    )
    if (
        constants[Constants.SHORT_STEP_THRESHOLD] <= 0.0
        or constants[Constants.SHORT_STEP_THRESHOLD] >= 1.0
    ):
        raise ValueError(
            "The constant short_step_threshold must be in the interval (0, 1)."
        )
    constants.setdefault(
        Constants.LOW_RADIUS_FACTOR.value,
        DEFAULT_CONSTANTS[Constants.LOW_RADIUS_FACTOR],
    )
    constants[Constants.LOW_RADIUS_FACTOR.value] = float(
        constants[Constants.LOW_RADIUS_FACTOR]
    )
    if (
        constants[Constants.LOW_RADIUS_FACTOR] <= 0.0
        or constants[Constants.LOW_RADIUS_FACTOR] >= 1.0
    ):
        raise ValueError(
            "The constant low_radius_factor must be in the interval (0, 1)."
        )
    constants.setdefault(
        Constants.BYRD_OMOJOKUN_FACTOR.value,
        DEFAULT_CONSTANTS[Constants.BYRD_OMOJOKUN_FACTOR],
    )
    constants[Constants.BYRD_OMOJOKUN_FACTOR.value] = float(
        constants[Constants.BYRD_OMOJOKUN_FACTOR]
    )
    if (
        constants[Constants.BYRD_OMOJOKUN_FACTOR] <= 0.0
        or constants[Constants.BYRD_OMOJOKUN_FACTOR] >= 1.0
    ):
        raise ValueError(
            "The constant byrd_omojokun_factor must be in the interval (0, 1)."
        )
    constants.setdefault(
        Constants.THRESHOLD_RATIO_CONSTRAINTS.value,
        DEFAULT_CONSTANTS[Constants.THRESHOLD_RATIO_CONSTRAINTS],
    )
    constants[Constants.THRESHOLD_RATIO_CONSTRAINTS.value] = float(
        constants[Constants.THRESHOLD_RATIO_CONSTRAINTS]
    )
    if constants[Constants.THRESHOLD_RATIO_CONSTRAINTS] <= 1.0:
        raise ValueError(
            "The constant threshold_ratio_constraints must be greater than 1."
        )
    constants.setdefault(
        Constants.LARGE_SHIFT_FACTOR.value,
        DEFAULT_CONSTANTS[Constants.LARGE_SHIFT_FACTOR],
    )
    constants[Constants.LARGE_SHIFT_FACTOR.value] = float(
        constants[Constants.LARGE_SHIFT_FACTOR]
    )
    if constants[Constants.LARGE_SHIFT_FACTOR] < 0.0:
        raise ValueError("The constant large_shift_factor must be "
                         "nonnegative.")
    constants.setdefault(
        Constants.LARGE_GRADIENT_FACTOR.value,
        DEFAULT_CONSTANTS[Constants.LARGE_GRADIENT_FACTOR],
    )
    constants[Constants.LARGE_GRADIENT_FACTOR.value] = float(
        constants[Constants.LARGE_GRADIENT_FACTOR]
    )
    if constants[Constants.LARGE_GRADIENT_FACTOR] <= 1.0:
        raise ValueError(
            "The constant large_gradient_factor must be greater than 1."
        )
    constants.setdefault(
        Constants.RESOLUTION_FACTOR.value,
        DEFAULT_CONSTANTS[Constants.RESOLUTION_FACTOR],
    )
    constants[Constants.RESOLUTION_FACTOR.value] = float(
        constants[Constants.RESOLUTION_FACTOR]
    )
    if constants[Constants.RESOLUTION_FACTOR] <= 1.0:
        raise ValueError(
            "The constant resolution_factor must be greater than 1."
        )
    constants.setdefault(
        Constants.IMPROVE_TCG.value,
        DEFAULT_CONSTANTS[Constants.IMPROVE_TCG],
    )
    constants[Constants.IMPROVE_TCG.value] = bool(
        constants[Constants.IMPROVE_TCG]
    )

    # Check whether they are any unknown options.
    for key in kwargs:
        if key not in Constants.__members__.values():
            warnings.warn(f"Unknown constant: {key}.", RuntimeWarning, 3)
    return constants


def _eval(pb, framework, step, options):
    """
    Evaluate the objective and constraint functions.
    """
    if pb.n_eval >= options[Options.MAX_EVAL]:
        raise MaxEvalError
    x_eval = framework.x_best + step
    fun_val, cub_val, ceq_val = pb(x_eval, framework.penalty)
    r_val = pb.maxcv(x_eval, cub_val, ceq_val)
    if (
        fun_val <= options[Options.TARGET]
        and r_val <= options[Options.FEASIBILITY_TOL]
    ):
        raise TargetSuccess
    if pb.is_feasibility and r_val <= options[Options.FEASIBILITY_TOL]:
        raise FeasibleSuccess
    return fun_val, cub_val, ceq_val


def _build_result(pb, penalty, success, status, n_iter, options):
    """
    Build the result of the optimization process.
    """
    # Build the result.
    x, fun, maxcv = pb.best_eval(penalty)
    success = success and np.isfinite(fun) and np.isfinite(maxcv)
    if status not in [ExitStatus.TARGET_SUCCESS, ExitStatus.FEASIBLE_SUCCESS]:
        success = success and maxcv <= options[Options.FEASIBILITY_TOL]
    result = OptimizeResult()
    result.message = {
        ExitStatus.RADIUS_SUCCESS: "The lower bound for the trust-region "
                                   "radius has been reached",
        ExitStatus.TARGET_SUCCESS: "The target objective function value has "
                                   "been reached",
        ExitStatus.FIXED_SUCCESS: "All variables are fixed by the bound "
                                  "constraints",
        ExitStatus.CALLBACK_SUCCESS: "The callback requested to stop the "
                                     "optimization procedure",
        ExitStatus.FEASIBLE_SUCCESS: "The feasibility problem received has "
                                     "been solved successfully",
        ExitStatus.MAX_EVAL_WARNING: "The maximum number of function "
                                     "evaluations has been exceeded",
        ExitStatus.MAX_ITER_WARNING: "The maximum number of iterations has "
                                     "been exceeded",
        ExitStatus.INFEASIBLE_ERROR: "The bound constraints are infeasible",
        ExitStatus.LINALG_ERROR: "A linear algebra error occurred",
    }.get(status, "Unknown exit status")
    result.success = success
    result.status = status.value
    result.x = pb.build_x(x)
    result.fun = fun
    result.maxcv = maxcv
    result.nfev = pb.n_eval
    result.nit = n_iter
    if options[Options.STORE_HISTORY]:
        result.fun_history = pb.fun_history
        result.maxcv_history = pb.maxcv_history

    # Print the result if requested.
    if options[Options.VERBOSE]:
        _print_step(
            result.message,
            pb,
            result.x,
            result.fun,
            result.maxcv,
            result.nfev,
            result.nit,
        )
    return result


def _print_step(message, pb, x, fun_val, r_val, n_eval, n_iter):
    """
    Print information about the current state of the optimization process.
    """
    print()
    print(f"{message}.")
    print(f"Number of function evaluations: {n_eval}.")
    print(f"Number of iterations: {n_iter}.")
    if not pb.is_feasibility:
        print(f"Least value of {pb.fun_name}: {fun_val}.")
    print(f"Maximum constraint violation: {r_val}.")
    with np.printoptions(**PRINT_OPTIONS):
        print(f"Corresponding point: {x}.")