# Bounds may appear unused in this file but we need to import it to make it available to the user from scipy.optimize import NonlinearConstraint, LinearConstraint, Bounds from .common._nonlinear_constraints import process_nl_constraints from .common._linear_constraints import ( combine_multiple_linear_constraints, separate_LC_into_eq_and_ineq, ) from .common._bounds import process_bounds from enum import Enum from .common._project import _project from .common.linalg import get_arrays_tol from .cobyla.cobyla import cobyla import numpy as np from collections.abc import Iterable class ConstraintType(Enum): LINEAR_OBJECT = 5 NONLINEAR_OBJECT = 10 LINEAR_DICT = 15 NONLINEAR_DICT = 20 def get_constraint_type(constraint): if isinstance(constraint, dict) and ("A" in constraint) and ("lb" in constraint) and ("ub" in constraint): return ConstraintType.LINEAR_DICT elif isinstance(constraint, dict) and ("fun" in constraint) and ("lb" in constraint) and ("ub" in constraint): return ConstraintType.NONLINEAR_DICT elif hasattr(constraint, "A") and hasattr(constraint, "lb") and hasattr(constraint, "ub"): return ConstraintType.LINEAR_OBJECT elif hasattr(constraint, "fun") and hasattr(constraint, "lb") and hasattr(constraint, "ub"): return ConstraintType.NONLINEAR_OBJECT else: raise ValueError(f"Constraint type {type(constraint)} not recognized") def process_constraints(constraints): # First throw it back if it's an empty tuple if not constraints: return None, None # Next figure out if it's a list of constraints or a single constraint # If it's a single constraint, make it a list, and then the remaining logic # doesn't have to change if not isinstance(constraints, Iterable): constraints = [constraints] # Separate out the linear and nonlinear constraints linear_constraints = [] nonlinear_constraints = [] for constraint in constraints: constraint_type = get_constraint_type(constraint) if constraint_type is ConstraintType.LINEAR_OBJECT: linear_constraints.append(constraint) elif constraint_type is ConstraintType.NONLINEAR_OBJECT: nonlinear_constraints.append(constraint) elif constraint_type == ConstraintType.LINEAR_DICT: linear_constraints.append(LinearConstraint(constraint["A"], constraint["lb"], constraint["ub"])) elif constraint_type == ConstraintType.NONLINEAR_DICT: nonlinear_constraints.append(NonlinearConstraint(constraint["fun"], constraint["lb"], constraint["ub"])) else: raise ValueError("Constraint type not recognized") if len(nonlinear_constraints) > 0: nonlinear_constraint_function = process_nl_constraints(nonlinear_constraints) else: nonlinear_constraint_function = None # Determine if we have multiple linear constraints, just 1, or none, and process accordingly if len(linear_constraints) > 1: linear_constraint = combine_multiple_linear_constraints(linear_constraints) elif len(linear_constraints) == 1: linear_constraint = linear_constraints[0] else: linear_constraint = None return linear_constraint, nonlinear_constraint_function def minimize(fun, x0, args=(), method=None, bounds=None, constraints=(), callback=None, options=None): linear_constraint, nonlinear_constraint_function = process_constraints(constraints) options = {'quiet': True} if options is None else options quiet = options.get("quiet", True) if method is None: if nonlinear_constraint_function is not None: if not quiet: print("Nonlinear constraints detected, applying COBYLA") method = "cobyla" elif linear_constraint is not None: if not quiet: print("Linear constraints detected without nonlinear constraints, applying LINCOA") method = "lincoa" elif bounds is not None: if not quiet: print("Bounds without linear or nonlinear constraints detected, applying BOBYQA") method = "bobyqa" else: if not quiet: print("No bounds or constraints detected, applying NEWUOA") method = "newuoa" else: # Raise some errors if methods were called with inappropriate options method = method.lower() if method not in ('newuoa', 'uobyqa', 'bobyqa', 'cobyla', 'lincoa'): raise ValueError(f"Method must be one of NEWUOA, UOBYQA, BOBYQA, COBYLA, or LINCOA, not '{method}'") if method != "cobyla" and nonlinear_constraint_function is not None: raise ValueError("Nonlinear constraints were provided for an algorithm that cannot handle them") if method not in ("cobyla", "lincoa") and linear_constraint is not None: raise ValueError("Linear constraints were provided for an algorithm that cannot handle them") if method not in ("cobyla", "bobyqa", "lincoa") and bounds is not None: raise ValueError("Bounds were provided for an algorithm that cannot handle them") # Try to get the length of x0. If we can't that likely means it's a scalar, and # in that case we turn it into an array and wrap the original function so that it # can accept an array and return a scalar. try: lenx0 = len(x0) except TypeError: x0 = np.array([x0]) original_scalar_fun = fun def scalar_fun(x): return original_scalar_fun(x[0], *args) fun = scalar_fun lenx0 = 1 lb, ub = process_bounds(bounds, lenx0) # Check which variables are fixed and eliminate them from the problem. # Save the indices and values so that we can call the original function with # an array of the appropriate size, and so that we can add the fixed values to the # result when COBYLA returns. tol = get_arrays_tol(lb, ub) _fixed_idx = ( (lb <= ub) & (np.abs(lb - ub) < tol) ) if any(_fixed_idx): _fixed_values = 0.5 * ( lb[_fixed_idx] + ub[_fixed_idx] ) _fixed_values = np.clip( _fixed_values, lb[_fixed_idx], ub[_fixed_idx], ) x0 = x0[~_fixed_idx] lb = lb[~_fixed_idx] ub = ub[~_fixed_idx] original_fun = fun def fixed_fun(x): newx = np.zeros(lenx0) newx[_fixed_idx] = _fixed_values newx[~_fixed_idx] = x return original_fun(newx, *args) fun = fixed_fun # Project x0 onto the feasible set if nonlinear_constraint_function is None: result = _project(x0, lb, ub, {"linear": linear_constraint, "nonlinear": None}) x0 = result.x if linear_constraint is not None: A_eq, b_eq, A_ineq, b_ineq = separate_LC_into_eq_and_ineq(linear_constraint) else: A_eq, b_eq, A_ineq, b_ineq = None, None, None, None if nonlinear_constraint_function is not None: # If there is a nonlinear constraint function, we will call COBYLA, which needs the number of nonlinear # constraints (m_nlcon). In order to get this number we need to evaluate the constraint function at x0. # The constraint value at x0 (nlconstr0) is not discarded but passed down to the Fortran backend, as its # evaluation is assumed to be expensive. We also evaluate the objective function at x0 and pass the result # (f0) down to the Fortran backend, which expects nlconstr0 and f0 to be provided in sync. def calcfc(x): f = fun(x, *args) nlconstr = nonlinear_constraint_function(x) return f, nlconstr else: def calcfc(x): f = fun(x, *args) constr = np.zeros(0) return f, constr f0, nlconstr0 = calcfc(x0) if 'quiet' in options: del options['quiet'] if 'maxfev' in options: options['maxfun'] = options['maxfev'] del options['maxfev'] result = cobyla( calcfc, len(nlconstr0), x0, A_ineq, b_ineq, A_eq, b_eq, lb, ub, f0=f0, nlconstr0=nlconstr0, callback=callback, **options ) if any(_fixed_idx): newx = np.zeros(lenx0) newx[_fixed_idx] = _fixed_values newx[~_fixed_idx] = result.x result.x = newx return result