diff --git a/CHANGES.rst b/CHANGES.rst index 20046df27..e8c4a7681 100644 --- a/CHANGES.rst +++ b/CHANGES.rst @@ -33,7 +33,7 @@ Breaking changes ``estimagic.OptimizeLogReader``. - Convenience functions to create namedtuples are removed from ``estimagic.utilities``. - +- :gh:`346` Add option to use nonlinear constraints (:ghuser:`timmens`) - :gh:`345` Moves estimation_table to new latex functionality of pandas (:ghuser:`mpetrosian`) - :gh:`344` Adds pytree support to slice_plot (:ghuser:`janosg`) diff --git a/docs/source/explanations/optimization/implementation_of_constraints.rst b/docs/source/explanations/optimization/implementation_of_constraints.rst index 5f236e3ba..2c017b04f 100644 --- a/docs/source/explanations/optimization/implementation_of_constraints.rst +++ b/docs/source/explanations/optimization/implementation_of_constraints.rst @@ -24,6 +24,14 @@ constraints into constrained optimizers. Reparametrization and penalties. Below explain what both approaches are, why we chose the reparametrization approach over penalties and which reparametrizations we are using for each type of constraint. +.. note:: + + In this text we focus on constraints that can solved by estimagic via bijective and + differentiable transformations. General nonlinear constraints do not fall into this + category. If you want to use nonlinear constraints you can still do so, but + estimagic will simply pass the constraints to your chosen optimizer. See + :ref:`constraints` for more details. + Possible approaches ------------------- @@ -32,9 +40,9 @@ Possible approaches Reparametrizations ~~~~~~~~~~~~~~~~~~ -In the reparametrization approach need to find an invertible mapping :math:`g` such as -well as two :math:`k'` dimensional vectors :math:`l` and :math:`u` such that: - +In the reparametrization approach we need to find an invertible mapping +:math:`g : \mathbb{R}^{k'} \to \mathbb{R}^k`, and two new bounds :math:`l'` and +:math:`u'` such that: .. math:: @@ -175,7 +183,7 @@ A suitable choice of :math:`\mathbf{\tilde{X}}` and :math:`\mathbf{M}` are: .. math:: \mathbf{\tilde{X}} \equiv \{(\tilde{x}_1, \tilde{x}_2)^T \mid \mathbf{\tilde{x}}_1 - \in \mathbb{R}^{k}$ \text{ and } \mathbf{l} \leq \mathbf{\tilde{x}}_2 \leq \mathbf{l}\} + \in \mathbb{R}^{k} \text{ and } \mathbf{l} \leq \mathbf{\tilde{x}}_2 \leq \mathbf{l}\} \mathbf{M} = \left[ {\begin{array}{cc} diff --git a/docs/source/explanations/optimization/internal_optimizers.rst b/docs/source/explanations/optimization/internal_optimizers.rst index 9b036588d..933b326ad 100644 --- a/docs/source/explanations/optimization/internal_optimizers.rst +++ b/docs/source/explanations/optimization/internal_optimizers.rst @@ -43,6 +43,7 @@ this interface. The mandatory conditions for an internal optimizer function are: first derivative jointly - lower_bounds: for lower bounds in form of a 1d numpy array - upper_bounds: for upper bounds in form of a 1d numpy array + - nonlinear_constraints: for nonlinear constraints in form a list of dictionaries Of course, algorithms that do not need a certain argument (e.g. unbounded or derivative free ones) do not need those arguments at all. @@ -98,6 +99,53 @@ different optimizers. However, complete transparency is possible and we try to d the exact meaning of all options for all optimizers. +Nonlinear constraints +--------------------- + +Estimagic can pass nonlinear constraints to the internal optimizer. The internal +interface for nonlinear constraints is as follows. + +A nonlinear constraint is a ``list`` of ``dict`` 's, where each ``dict`` represents a +group of constraints. In each group the constraint function can potentially be +multi-dimensional. We distinguish between equality and inequality constraints, which is +signalled by a dict entry ``type`` that takes values ``"eq"`` and ``"ineq"``. The +constraint function, which takes as input an internal parameter vector, is stored under +the entry ``fun``, while the Jacobian of that function is stored at ``jac``. The +tolerance for the constraints is stored under ``tol``. At last, the number of +constraints in each group is specified under ``n_constr``. An example list with one +constraint that would be passed to the internal optimizer is given by + +.. code-block:: + + constraints = [ + { + "type": "ineq", + "n_constr": 1, + "tol": 1e-5, + "fun": lambda x: x**3, + "jac": lambda x: 3 * x**2, + } + ] + + +.. note:: + + **Equality.** Internal equality constraints assume that the constraint is met when the function is + zero. That is + + .. math:: + + 0 = g(x) \in \mathbb{R}^m . + + **Inequality.** Internal inequality constraints assume that the constraint is met when the function is + greater or equal to zero. That is + + .. math:: + + 0 \leq g(x) \in \mathbb{R}^m . + + + Other conventions ----------------- diff --git a/docs/source/how_to_guides/optimization/how_to_specify_constraints.rst b/docs/source/how_to_guides/optimization/how_to_specify_constraints.rst index 2c45a7c98..6281aeff5 100644 --- a/docs/source/how_to_guides/optimization/how_to_specify_constraints.rst +++ b/docs/source/how_to_guides/optimization/how_to_specify_constraints.rst @@ -10,25 +10,25 @@ Constraints vs bounds Estimagic distinguishes between bounds and constraints. Bounds are lower and upper bounds for parameters. In the literature they are sometimes called box constraints. -Bounds are specified as "lower_bound" and "upper_bound" column of a params DataFrame -or as pytrees via the ``lower_bounds`` and ``upper_bounds`` argument to ``maximize`` and -``minimize``. +Bounds are specified as ``lower_bounds`` and ``upper_bounds`` argument to ``maximize`` +and ``minimize``. Examples with bounds can be found in `this tutorial`_. .. _this tutorial: ../../getting_started/first_optimization_with_estimagic.ipynb -Constraints are more general constraints on the parameters. This ranges from rather -simple ones (e.g. Parameters are fixed to a value, a group of parameters is required -to be equal) to more complex ones (like general linear constraints). +To specify more general constraints on the parameters you use can use the argument +``constraints``. This ranges from rather simple ones (e.g. parameters are fixed to a +value, a group of parameters is required to be equal) to more complex ones (like general +linear constraints, or even nonlinear constraints). Can you use constraints with all optimizers? ============================================ -We implement constraints via reparametrizations. Details are explained `here`_. This -means that you can use all of the constraints with any optimizer that supports -bounds. Some constraints (e.g. fixing parameters) can even be used with optimizers -that do not support bounds. +With the exception of general nonlinear constraints, we implement constraints via +reparametrizations. Details are explained `here`_. This means that you can use all of +the constraints with any optimizer that supports bounds. Some constraints (e.g. fixing +parameters) can even be used with optimizers that do not support bounds. .. _here: ../../explanations/optimization/implementation_of_constraints.rst @@ -288,8 +288,8 @@ flat numpy array are explained in the next section. typically it is more convenient to use the special cases instead of expressing them as a linear constraint. Internally, it will make no difference. - Let's impose the constraint that the sum of the average of the first four parameters - is at least 0.95. + Let's impose the constraint that the average of the first four parameters is at + least 0.95. .. code-block:: python @@ -315,6 +315,44 @@ flat numpy array are explained in the next section. can also be arrays (or even pytrees) with bounds and weights for each selected parameter. +.. tabbed:: nonlinear + + .. warning:: + + General nonlinear constraints that are specified via a black-box constraint + function can only be used if you choose an optimizer that supports it. + The feature is currently supported by the algorithms: + + * ``ipopt`` + * ``nlopt``: ``cobyla``, ``slsqp``, ``isres``, ``mma`` + * ``scipy``: ``cobyla``, ``slsqp``, ``trust_constr`` + + You can use nonlinear constraints to express restrictions of the form + ``lower_bound <= func(x) <= upper_bound`` or + ``func(x) = value`` where ``x`` are the selected parameters and ``func`` is the + constraint function. + + Let's impose the constraint that the product of all but the last parameter is 1. + + .. code-block:: python + + res = minimize( + criterion=criterion, + params=np.ones(6), + algorithm="scipy_slsqp", + constraints={ + "type": "nonlinear", + "selector": lambda x: x[:-1], + "func": lambda x: np.prod(x), + "value": 1.0, + }, + ) + + This yields: + + >>> array([ 1.31, 1.16, 1.01, 0.87, 0.75, -0. ]) + + Where the product of the all but the last parameters is equal to 1. Imposing multiple constraints at once diff --git a/src/estimagic/optimization/algo_options.py b/src/estimagic/optimization/algo_options.py index c38639d25..9a2b0c48f 100644 --- a/src/estimagic/optimization/algo_options.py +++ b/src/estimagic/optimization/algo_options.py @@ -179,6 +179,13 @@ """ +CONSTRAINTS_ABSOLUTE_TOLERANCE = 1e-5 +"""float: Allowed tolerance of the equality and inequality constraints for values to be +considered 'feasible'. + +""" + + """ ------------------------- Trust Region Parameters diff --git a/src/estimagic/optimization/cyipopt_optimizers.py b/src/estimagic/optimization/cyipopt_optimizers.py index 125df0217..49a91ed98 100644 --- a/src/estimagic/optimization/cyipopt_optimizers.py +++ b/src/estimagic/optimization/cyipopt_optimizers.py @@ -24,6 +24,8 @@ def ipopt( lower_bounds, upper_bounds, *, + # nonlinear constraints + nonlinear_constraints=(), # convergence criteria convergence_relative_criterion_tolerance=CONVERGENCE_RELATIVE_CRITERION_TOLERANCE, dual_inf_tol=1.0, @@ -494,7 +496,7 @@ def ipopt( x0=x, bounds=_get_scipy_bounds(lower_bounds, upper_bounds), jac=derivative, - constraints=(), + constraints=nonlinear_constraints, tol=convergence_relative_criterion_tolerance, options=options, ) diff --git a/src/estimagic/optimization/get_algorithm.py b/src/estimagic/optimization/get_algorithm.py index 1add6d913..1b853bee6 100644 --- a/src/estimagic/optimization/get_algorithm.py +++ b/src/estimagic/optimization/get_algorithm.py @@ -49,6 +49,7 @@ def get_final_algorithm( valid_kwargs, lower_bounds, upper_bounds, + nonlinear_constraints, algo_options, logging, db_kwargs, @@ -64,6 +65,8 @@ def get_final_algorithm( algorithm function. lower_bounds (np.ndarray): 1d numpy array with lower bounds. upper_bounds (np.ndarray): 1d numpy array with upper bounds. + nonlinear_constraints (list[dict]): List of dictionaries, each containing the + specification of a nonlinear constraint. algo_options (dict): Dictionary with additional keyword arguments for the algorithm. Entries that are not used by the algorithm are ignored with a warning. @@ -80,6 +83,7 @@ def get_final_algorithm( algo_options=algo_options, lower_bounds=lower_bounds, upper_bounds=upper_bounds, + nonlinear_constraints=nonlinear_constraints, algo_name=algo_name, valid_kwargs=valid_kwargs, ) @@ -253,6 +257,7 @@ def _adjust_options_to_algorithm( algo_options, lower_bounds, upper_bounds, + nonlinear_constraints, algo_name, valid_kwargs, ): @@ -291,4 +296,7 @@ def _adjust_options_to_algorithm( if "upper_bounds" in valid_kwargs: reduced["upper_bounds"] = upper_bounds + if "nonlinear_constraints" in valid_kwargs: + reduced["nonlinear_constraints"] = nonlinear_constraints + return reduced diff --git a/src/estimagic/optimization/nlopt_optimizers.py b/src/estimagic/optimization/nlopt_optimizers.py index 73443fc8a..e1af65fde 100644 --- a/src/estimagic/optimization/nlopt_optimizers.py +++ b/src/estimagic/optimization/nlopt_optimizers.py @@ -14,6 +14,9 @@ from estimagic.optimization.algo_options import ( STOPPING_MAX_CRITERION_EVALUATIONS_GLOBAL, ) +from estimagic.parameters.nonlinear_constraints import ( + equality_as_inequality_constraints, +) if IS_NLOPT_INSTALLED: import nlopt @@ -151,6 +154,7 @@ def nlopt_cobyla( lower_bounds, upper_bounds, *, + nonlinear_constraints=(), convergence_relative_params_tolerance=CONVERGENCE_RELATIVE_PARAMS_TOLERANCE, convergence_absolute_params_tolerance=CONVERGENCE_ABSOLUTE_PARAMS_TOLERANCE, convergence_relative_criterion_tolerance=CONVERGENCE_RELATIVE_CRITERION_TOLERANCE, @@ -170,6 +174,7 @@ def nlopt_cobyla( upper_bounds, algorithm=nlopt.LN_COBYLA, derivative=None, + nonlinear_constraints=nonlinear_constraints, convergence_xtol_rel=convergence_relative_params_tolerance, convergence_xtol_abs=convergence_absolute_params_tolerance, convergence_ftol_rel=convergence_relative_criterion_tolerance, @@ -403,6 +408,7 @@ def nlopt_mma( lower_bounds, upper_bounds, *, + nonlinear_constraints=(), convergence_relative_params_tolerance=CONVERGENCE_RELATIVE_PARAMS_TOLERANCE, convergence_absolute_params_tolerance=CONVERGENCE_ABSOLUTE_PARAMS_TOLERANCE, convergence_relative_criterion_tolerance=CONVERGENCE_RELATIVE_CRITERION_TOLERANCE, @@ -415,6 +421,9 @@ def nlopt_mma( For details see :ref:`list_of_nlopt_algorithms`. """ + # cannot handle equality constraints + nonlinear_constraints = equality_as_inequality_constraints(nonlinear_constraints) + out = _minimize_nlopt( criterion, x, @@ -422,6 +431,7 @@ def nlopt_mma( upper_bounds, algorithm=nlopt.LD_MMA, derivative=derivative, + nonlinear_constraints=nonlinear_constraints, convergence_xtol_rel=convergence_relative_params_tolerance, convergence_xtol_abs=convergence_absolute_params_tolerance, convergence_ftol_rel=convergence_relative_criterion_tolerance, @@ -492,6 +502,7 @@ def nlopt_slsqp( lower_bounds, upper_bounds, *, + nonlinear_constraints=(), convergence_relative_params_tolerance=CONVERGENCE_RELATIVE_PARAMS_TOLERANCE, convergence_absolute_params_tolerance=CONVERGENCE_ABSOLUTE_PARAMS_TOLERANCE, convergence_relative_criterion_tolerance=CONVERGENCE_RELATIVE_CRITERION_TOLERANCE, @@ -510,6 +521,7 @@ def nlopt_slsqp( upper_bounds, algorithm=nlopt.LD_SLSQP, derivative=derivative, + nonlinear_constraints=nonlinear_constraints, convergence_xtol_rel=convergence_relative_params_tolerance, convergence_xtol_abs=convergence_absolute_params_tolerance, convergence_ftol_rel=convergence_relative_criterion_tolerance, @@ -631,6 +643,7 @@ def nlopt_isres( lower_bounds, upper_bounds, *, + nonlinear_constraints=(), convergence_relative_params_tolerance=CONVERGENCE_RELATIVE_PARAMS_TOLERANCE, convergence_absolute_params_tolerance=CONVERGENCE_ABSOLUTE_PARAMS_TOLERANCE, convergence_relative_criterion_tolerance=CONVERGENCE_RELATIVE_CRITERION_TOLERANCE, @@ -648,6 +661,7 @@ def nlopt_isres( lower_bounds, upper_bounds, algorithm=nlopt.GN_ISRES, + nonlinear_constraints=nonlinear_constraints, convergence_xtol_rel=convergence_relative_params_tolerance, convergence_xtol_abs=convergence_absolute_params_tolerance, convergence_ftol_rel=convergence_relative_criterion_tolerance, @@ -714,6 +728,7 @@ def _minimize_nlopt( algorithm, *, derivative=None, + nonlinear_constraints=(), convergence_xtol_rel=None, convergence_xtol_abs=None, convergence_ftol_rel=None, @@ -748,11 +763,46 @@ def func(x, grad): opt.set_maxeval(stopping_max_eval) if population_size is not None: opt.set_population(population_size) + for constr in _get_nlopt_constraints(nonlinear_constraints, filter_type="eq"): + opt.add_equality_mconstraint(constr["fun"], constr["tol"]) + for constr in _get_nlopt_constraints(nonlinear_constraints, filter_type="ineq"): + opt.add_inequality_mconstraint(constr["fun"], constr["tol"]) opt.set_min_objective(func) solution_x = opt.optimize(x) return _process_nlopt_results(opt, solution_x) +def _get_nlopt_constraints(constraints, filter_type): + """Transform internal nonlinear constraints to NLOPT readable format.""" + filtered = [c for c in constraints if c["type"] == filter_type] + nlopt_constraints = [_internal_to_nlopt_constaint(c) for c in filtered] + return nlopt_constraints + + +def _internal_to_nlopt_constaint(c): + """ + Sign flip description: + + In estimagic, inequality constraints are internally defined as g(x) >= 0. NLOPT uses + h(x) <= 0, which is why we need to flip the sign. + + """ + tol = c["tol"] + if np.isscalar(tol): + tol = np.tile(tol, c["n_constr"]) + + def _constraint(result, x, grad): + result[:] = -c["fun"](x) # see docstring for sign flip + if grad.size > 0: + grad[:] = -c["jac"](x) # see docstring for sign flip + + new_constr = { + "fun": _constraint, + "tol": tol, + } + return new_constr + + def _process_nlopt_results(nlopt_obj, solution_x): messages = { 1: "Convergence achieved ", diff --git a/src/estimagic/optimization/optimize.py b/src/estimagic/optimization/optimize.py index 5867e6c3c..4d6a28864 100644 --- a/src/estimagic/optimization/optimize.py +++ b/src/estimagic/optimization/optimize.py @@ -24,6 +24,7 @@ from estimagic.optimization.tiktak import WEIGHT_FUNCTIONS from estimagic.parameters.conversion import aggregate_func_output_to_value from estimagic.parameters.conversion import get_converter +from estimagic.parameters.nonlinear_constraints import process_nonlinear_constraints from estimagic.process_user_function import process_func_of_params @@ -514,6 +515,22 @@ def _optimize( ) raise ValueError(msg.format(algo_info.name)) + # ================================================================================== + # Split constraints into nonlinear and reparametrization parts + # ================================================================================== + if isinstance(constraints, dict): + constraints = [constraints] + + nonlinear_constraints = [c for c in constraints if c["type"] == "nonlinear"] + + if nonlinear_constraints and "nonlinear_constraints" not in algo_kwargs: + raise ValueError( + f"Algorithm {algo_info.name} does not support nonlinear constraints." + ) + + # the following constraints will be handled via reparametrization + constraints = [c for c in constraints if c["type"] != "nonlinear"] + # ================================================================================== # prepare logging # ================================================================================== @@ -622,6 +639,7 @@ def _optimize( soft_upper_bounds=soft_upper_bounds, add_soft_bounds=multistart, ) + # ================================================================================== # initialize the log database # ================================================================================== @@ -662,6 +680,15 @@ def _optimize( direction=direction, ) + # process nonlinear constraints: + internal_constraints = process_nonlinear_constraints( + nonlinear_constraints=nonlinear_constraints, + params=params, + converter=converter, + numdiff_options=numdiff_options, + skip_checks=skip_checks, + ) + x = internal_params.values # ================================================================================== # get the internal algorithm @@ -672,6 +699,7 @@ def _optimize( valid_kwargs=algo_kwargs, lower_bounds=internal_params.lower_bounds, upper_bounds=internal_params.upper_bounds, + nonlinear_constraints=internal_constraints, algo_options=algo_options, logging=logging, db_kwargs=db_kwargs, diff --git a/src/estimagic/optimization/scipy_optimizers.py b/src/estimagic/optimization/scipy_optimizers.py index eafe26635..ab461d23f 100644 --- a/src/estimagic/optimization/scipy_optimizers.py +++ b/src/estimagic/optimization/scipy_optimizers.py @@ -59,8 +59,12 @@ from estimagic.optimization.algo_options import MAX_LINE_SEARCH_STEPS from estimagic.optimization.algo_options import STOPPING_MAX_CRITERION_EVALUATIONS from estimagic.optimization.algo_options import STOPPING_MAX_ITERATIONS +from estimagic.parameters.nonlinear_constraints import ( + equality_as_inequality_constraints, +) from estimagic.utilities import calculate_trustregion_initial_radius from scipy.optimize import Bounds +from scipy.optimize import NonlinearConstraint @mark_minimizer(name="scipy_lbfgsb") @@ -102,7 +106,7 @@ def scipy_lbfgsb( return process_scipy_result(res) -@mark_minimizer(name="scipy_lbfgsb") +@mark_minimizer(name="scipy_slsqp") def scipy_slsqp( criterion, derivative, @@ -110,6 +114,7 @@ def scipy_slsqp( lower_bounds, upper_bounds, *, + nonlinear_constraints=(), convergence_absolute_criterion_tolerance=CONVERGENCE_SECOND_BEST_ABSOLUTE_CRITERION_TOLERANCE, # noqa: E501 stopping_max_iterations=STOPPING_MAX_ITERATIONS, ): @@ -131,6 +136,7 @@ def scipy_slsqp( method="SLSQP", jac=derivative, bounds=_get_scipy_bounds(lower_bounds, upper_bounds), + constraints=nonlinear_constraints, options=options, ) @@ -304,6 +310,7 @@ def scipy_cobyla( criterion, x, *, + nonlinear_constraints=(), stopping_max_iterations=STOPPING_MAX_ITERATIONS, convergence_relative_params_tolerance=CONVERGENCE_RELATIVE_PARAMS_TOLERANCE, trustregion_initial_radius=None, @@ -318,10 +325,14 @@ def scipy_cobyla( options = {"maxiter": stopping_max_iterations, "rhobeg": trustregion_initial_radius} + # cannot handle equality constraints + nonlinear_constraints = equality_as_inequality_constraints(nonlinear_constraints) + res = scipy.optimize.minimize( fun=criterion, x0=x, method="COBYLA", + constraints=nonlinear_constraints, options=options, tol=convergence_relative_params_tolerance, ) @@ -390,6 +401,7 @@ def scipy_trust_constr( lower_bounds, upper_bounds, *, + nonlinear_constraints=(), convergence_absolute_gradient_tolerance=1e-08, convergence_relative_params_tolerance=CONVERGENCE_RELATIVE_PARAMS_TOLERANCE, stopping_max_iterations=STOPPING_MAX_ITERATIONS, @@ -410,12 +422,16 @@ def scipy_trust_constr( "initial_tr_radius": trustregion_initial_radius, } + # cannot handle equality constraints + nonlinear_constraints = equality_as_inequality_constraints(nonlinear_constraints) + res = scipy.optimize.minimize( fun=criterion_and_derivative, jac=True, x0=x, method="trust-constr", bounds=_get_scipy_bounds(lower_bounds, upper_bounds), + constraints=_get_scipy_constraints(nonlinear_constraints), options=options, ) @@ -444,6 +460,26 @@ def _get_scipy_bounds(lower_bounds, upper_bounds): return Bounds(lb=lower_bounds, ub=upper_bounds) +def _get_scipy_constraints(constraints): + """Transform internal nonlinear constraints to scipy readable format. + + This format is currently only used by scipy_trust-constr. + + """ + scipy_constraints = [_internal_to_scipy_constraint(c) for c in constraints] + return scipy_constraints + + +def _internal_to_scipy_constraint(c): + new_constr = NonlinearConstraint( + fun=c["fun"], + lb=np.zeros(c["n_constr"]), + ub=np.tile(np.inf, c["n_constr"]), + jac=c["jac"], + ) + return new_constr + + def _scipy_least_squares( criterion, derivative, diff --git a/src/estimagic/parameters/conversion.py b/src/estimagic/parameters/conversion.py index a6150bf81..108262784 100644 --- a/src/estimagic/parameters/conversion.py +++ b/src/estimagic/parameters/conversion.py @@ -143,8 +143,11 @@ def _params_from_internal(x, return_type="tree"): raise ValueError(msg) return out - def _derivative_to_internal(derivative_eval, x): - jacobian = tree_converter.derivative_flatten(derivative_eval) + def _derivative_to_internal(derivative_eval, x, jac_is_flat=False): + if jac_is_flat: + jacobian = derivative_eval + else: + jacobian = tree_converter.derivative_flatten(derivative_eval) x_unscaled = scale_converter.params_from_internal(x) jac_with_space_conversion = space_converter.derivative_to_internal( jacobian, x_unscaled @@ -224,10 +227,14 @@ def _fast_params_from_internal(x, return_type="tree"): def _get_fast_path_converter(params, lower_bounds, upper_bounds, primary_key): + def _fast_derivative_to_internal(derivative_eval, x, jac_is_flat=True): + # make signature compatible with non-fast path + return derivative_eval + converter = Converter( params_to_internal=lambda params: params.astype(float), params_from_internal=_fast_params_from_internal, - derivative_to_internal=lambda derivative_eval, x: derivative_eval, + derivative_to_internal=_fast_derivative_to_internal, func_to_internal=UNPACK_FUNCTIONS[primary_key], has_transforming_constraints=False, ) diff --git a/src/estimagic/parameters/nonlinear_constraints.py b/src/estimagic/parameters/nonlinear_constraints.py new file mode 100644 index 000000000..f2d09a139 --- /dev/null +++ b/src/estimagic/parameters/nonlinear_constraints.py @@ -0,0 +1,461 @@ +from functools import partial + +import numpy as np +import pandas as pd +from estimagic.differentiation.derivatives import first_derivative +from estimagic.exceptions import InvalidConstraintError +from estimagic.exceptions import InvalidFunctionError +from estimagic.optimization.algo_options import CONSTRAINTS_ABSOLUTE_TOLERANCE +from estimagic.parameters.block_trees import block_tree_to_matrix +from estimagic.parameters.tree_registry import get_registry +from pybaum import tree_flatten +from pybaum import tree_just_flatten +from pybaum import tree_unflatten + + +def process_nonlinear_constraints( + nonlinear_constraints, + params, + converter, + numdiff_options, + skip_checks, +): + """Process and prepare nonlinear constraints for internal use. + + Args: + nonlinear_constraints (list[dict]): List of dictionaries, each representing a + nonlinear constraint. + params (pandas): A pytree containing the parameters with respect to which the + criterion is optimized. Examples are a numpy array, a pandas Series, + a DataFrame with "value" column, a float and any kind of (nested) dictionary + or list containing these elements. See :ref:`params` for examples. + converter (Converter): NamedTuple with methods to convert between internal and + external parameters, derivatives and function outputs. + numdiff_options (dict): Keyword arguments for the calculation of numerical + derivatives. See :ref:`first_derivative` for details. Note that the default + method is changed to "forward" for speed reasons. Contains lower and upper + bounds of parameters. + skip_checks (bool): Whether checks on the inputs are skipped. This makes the + optimization faster, especially for very fast constraint functions. Default + False. + + Returns: + list[dict]: List of processed constraints. + + """ + # do checks first to fail fast + constraint_evals = [] + for _constraint in nonlinear_constraints: + + _eval = _check_validity_and_return_evaluation(_constraint, params, skip_checks) + constraint_evals.append(_eval) + + processed = [] + for _constraint, _eval in zip(nonlinear_constraints, constraint_evals): + + _processed_constraint = _process_nonlinear_constraint( + _constraint, + constraint_eval=_eval, + params=params, + converter=converter, + numdiff_options=numdiff_options, + ) + processed.append(_processed_constraint) + + return processed + + +def _process_nonlinear_constraint( + c, constraint_eval, params, converter, numdiff_options +): + """Process a single nonlinear constraint.""" + + # ================================================================================== + # Process selector and evaluate functions if necessary + # ================================================================================== + + external_selector = _process_selector(c) # functional selector + + constraint_func = c["func"] + + if constraint_eval is None: + selected = external_selector(params) + constraint_eval = constraint_func(selected) + + _n_constr = len(np.atleast_1d(constraint_eval)) + + # ================================================================================== + # Consolidate and transform jacobian + # ================================================================================== + + # process numdiff_options for numerical derivative + options = numdiff_options.copy() + options.pop("lower_bounds", None) + options.pop("upper_bounds", None) + + if "derivative" in c: + if not callable(c["derivative"]): + msg = "Jacobian of constraints needs to be callable." + raise ValueError(msg) + jacobian = c["derivative"] + else: + # use finite-differences if no closed-form jacobian is defined + def jacobian(p): + return first_derivative( + constraint_func, + p, + **options, + )["derivative"] + + # To define the internal Jacobian we need to know which parameters enter the + # contraint function. + selection_indices, n_params = _get_selection_indices(params, external_selector) + + def _internal_jacobian(x): + """Return Jacobian of constraint at internal parameters. + + The constraint function is written to be evaluated on a selection of the + external parameters. The optimizer, however, only works on internal parameters. + These can be significantly different from the external parameters, due to + estimagic's reparametrization features. In this function we compute the Jacobian + of the constraint at the internal parameters using information on the Jacobian + of the constraint at the selected external parameters. + + """ + params = converter.params_from_internal(x) + selected = external_selector(params) + jac = jacobian(selected) + jac_matrix = block_tree_to_matrix(jac, constraint_eval, selected) + jac_extended = _extend_jacobian(jac_matrix, selection_indices, n_params) + jac_internal = converter.derivative_to_internal( + jac_extended, x, jac_is_flat=True + ) + return np.atleast_2d(jac_internal) + + # ================================================================================== + # Transform constraint function and derive bounds + # ================================================================================== + _type = "eq" if "value" in c else "ineq" + + if _type == "eq": + + # ============================================================================== + # Equality constraints + # + # We define the internal constraint function to be satisfied if it is equal + # to zero, by subtracting the fixed value. + + _value = np.atleast_1d(np.array(c["value"], dtype=float)) + + def internal_constraint_func(x): + params = converter.params_from_internal(x) + select = external_selector(params) + out = np.atleast_1d(constraint_func(select)) - _value + return out + + jacobian_from_internal = _internal_jacobian + n_constr = _n_constr + + else: + + # ============================================================================== + # Inequality constraints + # + # We define the internal constraint function to be satisfied if it is + # greater than or equal to zero (positivity constraint). If the bounds already + # satify this condition we do not change anything, otherwise we need to perform + # a transformation. + + def _internal_constraint_func(x): + params = converter.params_from_internal(x) + select = external_selector(params) + return np.atleast_1d(constraint_func(select)) + + lower_bounds = c.get("lower_bounds", 0) + upper_bounds = c.get("upper_bounds", np.inf) + + transformation = _get_transformation(lower_bounds, upper_bounds) + + internal_constraint_func = _compose_funcs( + _internal_constraint_func, transformation["func"] + ) + + jacobian_from_internal = _compose_funcs( + _internal_jacobian, transformation["derivative"] + ) + + n_constr = 2 * _n_constr if transformation["name"] == "stack" else _n_constr + + internal_constr = { + "n_constr": n_constr, + "type": _type, + "fun": internal_constraint_func, # internal name for 'func' + "jac": jacobian_from_internal, # internal name for 'derivative' + "tol": c.get("tol", CONSTRAINTS_ABSOLUTE_TOLERANCE), + } + + return internal_constr + + +def equality_as_inequality_constraints(nonlinear_constraints): + """Return constraints where equality constraints are converted to inequality.""" + constraints = [_equality_to_inequality(c) for c in nonlinear_constraints] + return constraints + + +def _equality_to_inequality(c): + """Transform a single constraint. + + An equality constaint g(x) = 0 can be transformed to two inequality constraints + using (g(x), -g(x)) >= 0. Hence, the number of constraints doubles, and the + constraint functions itself as well as the derivative need to be updated. + + """ + if c["type"] == "eq": + + def transform(x, func): + value = func(x) + return np.concatenate((value, -value), axis=0) + + out = { + "fun": partial(transform, func=c["fun"]), + "jac": partial(transform, func=c["jac"]), + "n_constr": 2 * c["n_constr"], + "tol": c["tol"], + "type": "ineq", + } + else: + out = c + return out + + +# ====================================================================================== +# Helper Functions +# ====================================================================================== + + +def _process_selector(c): + if "selector" in c: + selector = c["selector"] + elif "loc" in c: + + def selector(params): + return params.loc[c["loc"]] + + elif "query" in c: + + def selector(params): + return params.query(c["query"]) + + else: + selector = _identity + return selector + + +def _compose_funcs(f, g): + return lambda x: g(f(x)) + + +def _identity(x): + return x + + +# ====================================================================================== +# Jacobian helper functions +# ====================================================================================== + + +def _extend_jacobian(jac_mat, selection_indices, n_params): + """Extend Jacobian on selected parameters to full params. + + Jacobian of constraints is defined on a selection of the parameters, however, we + need the Jacobian on the full params. Since the Jacobian is trivially zero at the + non-selected params we can simply fill a zero matrix. + + """ + jac_extended = np.zeros((jac_mat.shape[0], n_params)) + jac_extended[:, selection_indices] = jac_mat + return jac_extended + + +def _get_selection_indices(params, selector): + """Get index of selected flat params and number of flat params.""" + registry = get_registry(extended=True) + flat_params, params_treedef = tree_flatten(params, registry=registry) + n_params = len(flat_params) + indices = np.arange(n_params, dtype=int) + params_indices = tree_unflatten(params_treedef, indices, registry=registry) + selected = selector(params_indices) + selection_indices = np.array( + tree_just_flatten(selected, registry=registry), dtype=int + ) + return selection_indices, n_params + + +# ====================================================================================== +# Transformation helper functions +# ====================================================================================== + + +def _get_transformation(lower_bounds, upper_bounds): + """Get transformation given bounds. + + The internal inequality constraint is defined as h(x) >= 0. However, the user can + specify: a <= g(x) <= b. To get the internal represenation we need to transform the + constraint. + + """ + transformation_type = _get_transformation_type(lower_bounds, upper_bounds) + + if transformation_type == "identity": + transformer = {"func": _identity, "derivative": _identity} + elif transformation_type == "subtract_lb": + transformer = { + "func": lambda v: v - lower_bounds, + "derivative": _identity, + } + elif transformation_type == "stack": + transformer = { + "func": lambda v: np.concatenate( + (v - lower_bounds, upper_bounds - v), axis=0 + ), + "derivative": lambda v: np.concatenate((v, -v), axis=0), + } + transformer["name"] = transformation_type + return transformer + + +def _get_transformation_type(lower_bounds, upper_bounds): + lb_is_zero = not np.count_nonzero(lower_bounds) + ub_is_inf = np.all(np.isposinf(upper_bounds)) + + if lb_is_zero and ub_is_inf: + # the external constraint is already in the correct format + _transformation_type = "identity" + elif ub_is_inf: + # the external constraint can be transformed by subtraction + _transformation_type = "subtract_lb" + else: + # the external constraint can only be transformed by duplication (stacking) + _transformation_type = "stack" + return _transformation_type + + +# ====================================================================================== +# Checks +# ====================================================================================== + + +def _check_validity_and_return_evaluation(c, params, skip_checks): + """Check that nonlinear constraints are valid. + + Returns: + constaint_eval: Evaluation of constraint at params, if skip_checks if False, + else None. + + """ + # ================================================================================== + # check functions + # ================================================================================== + + if "func" not in c: + raise InvalidConstraintError( + "Constraint needs to have entry 'fun', representing the constraint " + "function." + ) + if not callable(c["func"]): + raise InvalidConstraintError( + "Entry 'fun' in nonlinear constraints has be callable." + ) + + if "derivative" in c and not callable(c["derivative"]): + raise InvalidConstraintError( + "Entry 'jac' in nonlinear constraints has be callable." + ) + + # ================================================================================== + # check bounds + # ================================================================================== + + is_equality_constraint = "value" in c + + if is_equality_constraint: + if "lower_bounds" in c or "upper_bounds" in c: + raise InvalidConstraintError( + "Only one of 'value' or ('lower_bounds', 'upper_bounds') can be " + "passed to a nonlinear constraint." + ) + + if not is_equality_constraint: + if "lower_bounds" not in c and "upper_bounds" not in c: + raise InvalidConstraintError( + "For inequality constraint at least one of ('lower_bounds', " + "'upper_bounds') has to be passed to the nonlinear constraint." + ) + + if "lower_bounds" in c and "upper_bounds" in c: + if not np.all(np.array(c["lower_bounds"]) <= np.array(c["upper_bounds"])): + raise InvalidConstraintError( + "If lower bounds need to less than or equal to upper bounds." + ) + + # ================================================================================== + # check selector + # ================================================================================== + + if "selector" in c: + if not callable(c["selector"]): + raise InvalidConstraintError( + f"'selector' entry needs to be callable in constraint {c}." + ) + else: + try: + c["selector"](params) + except Exception: + raise InvalidFunctionError( + "Error when calling 'selector' function on params in constraint " + f" {c}" + ) + + elif "loc" in c: + if not isinstance(params, (pd.Series, pd.DataFrame)): + raise InvalidConstraintError( + "params needs to be pd.Series or pd.DataFrame to use 'loc' selector in " + f"in consrtaint {c}." + ) + try: + params.loc[c["loc"]] + except (KeyError, IndexError): + raise InvalidConstraintError("'loc' string is invalid.") + + elif "query" in c: + if not isinstance(params, pd.DataFrame): + raise InvalidConstraintError( + "params needs to be pd.DataFrame to use 'query' selector in " + f"constraints {c}." + ) + try: + params.query(c["query"]) + except Exception: + raise InvalidConstraintError( + f"'query' string is invalid in constraint {c}." + ) + + # ================================================================================== + # check that constraints can be evaluated + # ================================================================================== + + constraint_eval = None + + if not skip_checks: + + selector = _process_selector(c) + + try: + constraint_eval = c["func"](selector(params)) + except Exception: + raise InvalidFunctionError( + f"Error when evaluating function of constraint {c}." + ) + + return constraint_eval diff --git a/src/estimagic/parameters/process_constraints.py b/src/estimagic/parameters/process_constraints.py index 841d639b8..491558279 100644 --- a/src/estimagic/parameters/process_constraints.py +++ b/src/estimagic/parameters/process_constraints.py @@ -140,7 +140,7 @@ def _process_linear_weights(constraints): """Harmonize the weights of linear constraints. Args: - pc (list): Constraints where the selectors have already been processed. + constraints (list): Constraints where the selectors have already been processed. Returns: list: Constraints where all weights are Series. diff --git a/tests/optimization/test_with_nonlinear_constraints.py b/tests/optimization/test_with_nonlinear_constraints.py new file mode 100644 index 000000000..605f17729 --- /dev/null +++ b/tests/optimization/test_with_nonlinear_constraints.py @@ -0,0 +1,249 @@ +import itertools +import warnings + +import numpy as np +import pytest +from estimagic import maximize +from estimagic import minimize +from estimagic.optimization import AVAILABLE_ALGORITHMS +from numpy.testing import assert_array_almost_equal as aaae + + +NLC_ALGORITHMS = [ + name + for name, algo in AVAILABLE_ALGORITHMS.items() + if "nonlinear_constraints" in algo._algorithm_info.arguments +] + +# ====================================================================================== +# Two-dimension example with equality and inequality constraints +# ====================================================================================== + + +@pytest.fixture() +def nlc_2d_example(): + """Non-linear constraints: 2-dimensional example. + + See the example section in https://en.wikipedia.org/wiki/Nonlinear_programming. + + """ + + def criterion(x): + return np.sum(x) + + def derivative(x): + return np.ones_like(x) + + def constraint_func(x): + value = np.dot(x, x) + return np.array([value - 1, 2 - value]) + + def constraint_jac(x): + return 2 * np.row_stack((x.reshape(1, -1), -x.reshape(1, -1))) + + constraints_long = [ + { + "type": "nonlinear", + "func": constraint_func, + "derivative": constraint_jac, + "lower_bounds": np.zeros(2), + "tol": 1e-8, + } + ] + + constraints_flat = [ + { + "type": "nonlinear", + "func": lambda x: np.dot(x, x), + "derivative": lambda x: 2 * x, + "lower_bounds": 1, + "upper_bounds": 2, + "tol": 1e-8, + } + ] + + constraints_equality = [ + { + "type": "nonlinear", + "func": lambda x: np.dot(x, x), + "derivative": lambda x: 2 * x, + "value": 2, + } + ] + + constraints_equality_and_inequality = [ + { + "type": "nonlinear", + "func": lambda x: np.dot(x, x), + "derivative": lambda x: 2 * x, + "lower_bounds": 1, + }, + { + "type": "nonlinear", + "func": lambda x: np.dot(x, x), + "derivative": lambda x: 2 * x, + "value": 2, + }, + ] + + _kwargs = { + "criterion": criterion, + "params": np.array([0, np.sqrt(2)]), + "derivative": derivative, + "lower_bounds": np.zeros(2), + "upper_bounds": 2 * np.ones(2), + } + + kwargs = { + "flat": {**_kwargs, **{"constraints": constraints_flat}}, + "long": {**_kwargs, **{"constraints": constraints_long}}, + "equality": {**_kwargs, **{"constraints": constraints_equality}}, + "equality_and_inequality": { + **_kwargs, + **{"constraints": constraints_equality_and_inequality}, + }, + } + + solution_x = np.ones(2) + + return solution_x, kwargs + + +TEST_CASES = list( + itertools.product( + NLC_ALGORITHMS, ["flat", "long", "equality", "equality_and_inequality"] + ) +) + + +@pytest.mark.parametrize("algorithm, constr_type", TEST_CASES) +def test_nonlinear_optimization(nlc_2d_example, algorithm, constr_type): + """Test that available nonlinear optimizers solve a nonlinear constraints problem. + + We test for the cases of "equality", "inequality" and "equality_and_inequality" + constraints. + + """ + if "equality" in constr_type and algorithm == "nlopt_mma": + pytest.skip(reason="Very slow and low accuracy.") + + solution_x, kwargs = nlc_2d_example + if algorithm == "scipy_cobyla": + del kwargs[constr_type]["lower_bounds"] + del kwargs[constr_type]["upper_bounds"] + + with warnings.catch_warnings(): + warnings.simplefilter("ignore") + result = maximize(algorithm=algorithm, **kwargs[constr_type]) + + if AVAILABLE_ALGORITHMS[algorithm]._algorithm_info.is_global: + decimal = 0 + else: + decimal = 4 + + aaae(result.params, solution_x, decimal=decimal) + + +# ====================================================================================== +# Documentation example +# ====================================================================================== + + +def criterion(params): + offset = np.linspace(1, 0, len(params)) + x = params - offset + return x @ x + + +@pytest.mark.parametrize("algorithm", NLC_ALGORITHMS) +def test_documentation_example(algorithm): + if algorithm in ("nlopt_mma", "ipopt"): + pytest.skip(reason="Slow.") + + kwargs = { + "lower_bounds": np.zeros(6), + "upper_bounds": 2 * np.ones(6), + } + + if algorithm == "scipy_cobyla": + del kwargs["lower_bounds"] + del kwargs["upper_bounds"] + + minimize( + criterion=criterion, + params=np.ones(6), + algorithm=algorithm, + constraints={ + "type": "nonlinear", + "selector": lambda x: x[:-1], + "func": lambda x: np.prod(x), + "value": 1.0, + }, + **kwargs + ) + + +# ====================================================================================== +# Test: selection + reparametrization constraint + nonlinear constraint +# ====================================================================================== + + +@pytest.fixture() +def general_example(): + + params = {"a": np.array([0.1, 0.3, 0.4, 0.2]), "b": np.array([1.5, 2])} + + def criterion(params): + weights = np.array([0, 1, 2, 3]) + return params["a"] @ weights + params["b"].sum() + + def selector_probability_constraint(params): + return params["a"] + + def selector_nonlinear_constraint(params): + return {"probs": params["a"][:3][::-1], "unnecessary": params["b"]} + + def constraint(selected): + return selected["probs"] @ selected["probs"] + + constraints = [ + {"type": "probability", "selector": selector_probability_constraint}, + { + "type": "nonlinear", + "selector": selector_nonlinear_constraint, + "upper_bounds": 0.8, + "func": constraint, + "tol": 0.01, + }, + ] + + lower_bounds = {"b": np.array([0, 0])} + upper_bounds = {"b": np.array([2, 2])} + + kwargs = { + "criterion": criterion, + "params": params, + "constraints": constraints, + "lower_bounds": lower_bounds, + "upper_bounds": upper_bounds, + } + return kwargs + + +TEST_CASES = list(itertools.product(["ipopt"], [True, False])) + + +@pytest.mark.parametrize("algorithm, skip_checks", TEST_CASES) +def test_general_example(general_example, algorithm, skip_checks): + + kwargs = general_example + + with warnings.catch_warnings(): + warnings.simplefilter("ignore") + res = minimize(algorithm=algorithm, skip_checks=skip_checks, **kwargs) + + optimal_p1 = 0.5 + np.sqrt(3 / 20) # can be derived analytically + optimal_p2 = 1 - optimal_p1 + + aaae(res.params["a"], np.array([optimal_p1, optimal_p2, 0, 0]), decimal=4) + aaae(res.params["b"], np.array([0.0, 0]), decimal=5) diff --git a/tests/parameters/test_nonlinear_constraints.py b/tests/parameters/test_nonlinear_constraints.py new file mode 100644 index 000000000..b8d417881 --- /dev/null +++ b/tests/parameters/test_nonlinear_constraints.py @@ -0,0 +1,258 @@ +import itertools +from dataclasses import dataclass + +import numpy as np +import pandas as pd +import pytest +from estimagic.exceptions import InvalidConstraintError +from estimagic.parameters.nonlinear_constraints import ( + _check_validity_and_return_evaluation, +) +from estimagic.parameters.nonlinear_constraints import _get_selection_indices +from estimagic.parameters.nonlinear_constraints import _get_transformation +from estimagic.parameters.nonlinear_constraints import _get_transformation_type +from estimagic.parameters.nonlinear_constraints import _process_selector +from estimagic.parameters.nonlinear_constraints import ( + equality_as_inequality_constraints, +) +from estimagic.parameters.nonlinear_constraints import process_nonlinear_constraints +from estimagic.parameters.tree_registry import get_registry +from numpy.testing import assert_array_equal +from pandas.testing import assert_frame_equal +from pybaum import tree_just_flatten + + +@dataclass +class Converter: + def params_from_internal(self, x): + return x + + def params_to_internal(self, params): + registry = get_registry(extended=True) + return np.array(tree_just_flatten(params, registry=registry)) + + +# ====================================================================================== +# _get_transformation_type +# ====================================================================================== +TEST_CASES = [ + (0, np.inf, "identity"), # (lower_bounds, upper_bounds, expected) + (-1, 2, "stack"), + (np.zeros(3), np.ones(3), "stack"), + (np.zeros(3), np.tile(np.inf, 3), "identity"), + (np.array([1, 2]), np.tile(np.inf, 2), "subtract_lb"), +] + + +@pytest.mark.parametrize("lower_bounds, upper_bounds, expected", TEST_CASES) +def test_get_transformation_type(lower_bounds, upper_bounds, expected): + got = _get_transformation_type(lower_bounds, upper_bounds) + assert got == expected + + +# ====================================================================================== +# _get_transformation +# ====================================================================================== + +TEST_CASES = [ + # (lower_bounds, upper_bounds, case, expected) # noqa: E800 + (0, 0, "func", {"name": "stack", "out": np.array([1, -1])}), + (1, 1, "func", {"name": "stack", "out": np.array([0, 0])}), + (0, 0, "derivative", {"name": "stack", "out": np.array([1, -1])}), + (1, 1, "derivative", {"name": "stack", "out": np.array([1, -1])}), + (1, np.inf, "func", {"name": "subtract_lb", "out": np.array([0])}), + (0, np.inf, "derivative", {"name": "identity", "out": np.array([1])}), +] + + +@pytest.mark.parametrize("lower_bounds, upper_bounds, case, expected", TEST_CASES) +def test_get_positivity_transform(lower_bounds, upper_bounds, case, expected): + transform = _get_transformation(lower_bounds, upper_bounds) + + got = transform[case](np.array([1])) + assert np.all(got == expected["out"]) + assert transform["name"] == expected["name"] + + +# ====================================================================================== +# _get_selection_indices +# ====================================================================================== + + +def test_get_selection_indices(): + + params = {"a": [0, 1, 2], "b": [3, 4, 5]} + selector = lambda p: p["a"] + + expected = np.array([0, 1, 2], dtype=int) + got_index, got_n_params = _get_selection_indices(params, selector) + + assert got_n_params == 6 + assert_array_equal(got_index, expected) + + +# ====================================================================================== +# _process_selector +# ====================================================================================== +TEST_CASES = [ + ({"selector": lambda x: x**2}, 10, 100), # (constraint, params, expected) + ({"loc": "a"}, pd.Series([0, 1], index=["a", "b"]), 0), + ( + {"query": "a == 1"}, + pd.DataFrame([[1], [0]], columns=["a"]), + pd.DataFrame([[1]], columns=["a"]), + ), +] + + +@pytest.mark.parametrize("constraint, params, expected", TEST_CASES) +def test_process_selector(constraint, params, expected): + _selector = _process_selector(constraint) + got = _selector(params) + + if isinstance(got, pd.DataFrame): + assert_frame_equal(got, expected) + else: + assert got == expected + + +# ====================================================================================== +# _check_validity_nonlinear_constraint +# ====================================================================================== +TEST_CASES = [ + {}, # no fun + {"func": 10}, # non-callable fun + {"func": lambda x: x, "derivative": 10}, # non-callable jac + {"func": lambda x: x}, # no bounds at all + { + "func": lambda x: x, + "value": 1, + "lower_bounds": 1, + }, # cannot have value and bounds + { + "func": lambda x: x, + "value": 1, + "upper_bounds": 1, + }, # cannot have value and bounds + {"func": lambda x: x}, # needs to have at least one bound + {"func": lambda x: x, "lower_bounds": 1, "upper_bounds": 0}, + {"func": lambda x: x, "selector": 10}, + {"func": lambda x: x, "loc": 10}, + {"func": lambda x: x, "query": 10}, +] + +TEST_CASES = list( + itertools.product(TEST_CASES, [np.arange(3), pd.DataFrame({"a": [0, 1, 2]})]) +) + + +@pytest.mark.parametrize("constraint, params", TEST_CASES) +def test_check_validity_nonlinear_constraint(constraint, params): + with pytest.raises(InvalidConstraintError): + _check_validity_and_return_evaluation(constraint, params, skip_checks=False) + + +def test_check_validity_nonlinear_constraint_correct_example(): + constr = { + "func": lambda x: x, + "derivative": lambda x: np.ones_like(x), + "lower_bounds": np.arange(4), + "selector": lambda x: x[:1], + } + _check_validity_and_return_evaluation( + constr, params=np.arange(4), skip_checks=False + ) + + +# ====================================================================================== +# equality_as_inequality_constraints +# ====================================================================================== +TEST_CASES = [ + ( + [ + { + "type": "ineq", + "fun": lambda x: np.array([x]), + "jac": lambda x: np.array([[1]]), + "n_constr": 1, + } + ], # constraints + "same", # expected + ), + ( + [ + { + "type": "ineq", + "fun": lambda x: np.array([x]), + "jac": lambda x: np.array([[1]]), + "n_constr": 1, + } + ], # constraints + [ + { + "type": "eq", + "fun": lambda x: np.array([x, -x]).reshape(-1, 1), + "jac": lambda x: np.array([[1], [-1]]), + "n_constr": 1, + } + ], # expected + ), +] + + +@pytest.mark.parametrize("constraints, expected", TEST_CASES) +def test_equality_as_inequality_constraints(constraints, expected): + got = equality_as_inequality_constraints(constraints) + if expected == "same": + assert got == constraints + + for g, c in zip(got, constraints): + if c["type"] == "eq": + assert g["n_constr"] == 2 * c["n_constr"] + assert g["type"] == "ineq" + + +# ====================================================================================== +# process_nonlinear_constraints +# ====================================================================================== + + +def test_process_nonlinear_constraints(): + + nonlinear_constraints = [ + {"type": "nonlinear", "func": lambda x: np.dot(x, x), "value": 1}, + { + "type": "nonlinear", + "func": lambda x: x, + "lower_bounds": -1, + "upper_bounds": 2, + }, + ] + + params = np.array([1.0]) + + converter = Converter() + + numdiff_options = {"lower_bounds": params, "upper_bounds": params} + + got = process_nonlinear_constraints( + nonlinear_constraints, params, converter, numdiff_options, skip_checks=False + ) + + expected = [ + {"type": "eq", "fun": lambda x: np.dot(x, x) - 1.0, "n_constr": 1}, + { + "type": "ineq", + "fun": lambda x: np.concatenate((x + 1.0, 2.0 - x), axis=0), + "n_constr": 2, + }, + ] + + for g, e in zip(got, expected): + assert g["type"] == e["type"] + assert g["n_constr"] == e["n_constr"] + for x in [0.1, 0.2, 1.2, -2.0]: + x = np.array([x]) + assert_array_equal(g["fun"](x), e["fun"](x)) + assert "jac" in g + assert "tol" in g diff --git a/tests/visualization/test_derivative_plot.py b/tests/visualization/test_derivative_plot.py index 74c44b461..d97aff643 100644 --- a/tests/visualization/test_derivative_plot.py +++ b/tests/visualization/test_derivative_plot.py @@ -1,5 +1,3 @@ -import itertools - import numpy as np import pandas as pd import pytest @@ -71,16 +69,6 @@ def test__select_eval_with_lowest_and_highest_step(): assert_array_equal(got, expected) -def _powerset(iterable): - s = list(iterable) - pset = itertools.chain.from_iterable( - itertools.combinations(s, r) for r in range(len(s) + 1) - ) - pset = [e for e in pset if len(e) > 0] - pset = list(map(lambda x: x if len(x) > 1 else x[0], pset)) - return pset - - def f1(x): y1 = np.sin(x[0]) + np.cos(x[1]) + x[2] return y1