From cfcdf157fb29090e490ea4c92436f9d43a4c7dac Mon Sep 17 00:00:00 2001 From: Nickolai Belakovski Date: Mon, 29 Sep 2025 10:19:47 -0400 Subject: [PATCH 1/5] Adjust linear and nonlinear constraints upon eliminating fixed bounds When we eliminate fixed bounds, we need to update not only the objective function, but also the linear constraints and the nonlinear constraint function. Tests have been added to ensure the added code works, and the tests have been demonstrated to fail when the new code is removed. Ref: #249 --- pyprima/src/pyprima/__init__.py | 15 +++++++++++ pyprima/tests/test_bounds.py | 47 ++++++++++++++++++++++++++++++++- 2 files changed, 61 insertions(+), 1 deletion(-) diff --git a/pyprima/src/pyprima/__init__.py b/pyprima/src/pyprima/__init__.py index 577e1d2661..1c7f54160c 100644 --- a/pyprima/src/pyprima/__init__.py +++ b/pyprima/src/pyprima/__init__.py @@ -151,6 +151,21 @@ def fixed_fun(x): newx[~_fixed_idx] = x return original_fun(newx, *args) fun = fixed_fun + # Adjust linear_constraint + if linear_constraint: + new_lb = linear_constraint.lb - linear_constraint.A[:, _fixed_idx] @ _fixed_values + new_ub = linear_constraint.ub - linear_constraint.A[:, _fixed_idx] @ _fixed_values + new_A = linear_constraint.A[:, ~_fixed_idx] + linear_constraint = LinearConstraint(new_A, new_lb, new_ub) + # Adjust nonlinear constraints + if nonlinear_constraint_function: + original_nlc_fun = nonlinear_constraint_function + def fixed_nlc_fun(x): + newx = np.zeros(lenx0) + newx[_fixed_idx] = _fixed_values + newx[~_fixed_idx] = x + return original_nlc_fun(newx, *args) + nonlinear_constraint_function = fixed_nlc_fun # Project x0 onto the feasible set diff --git a/pyprima/tests/test_bounds.py b/pyprima/tests/test_bounds.py index 72d74db7b7..80fda6a6d7 100644 --- a/pyprima/tests/test_bounds.py +++ b/pyprima/tests/test_bounds.py @@ -1,4 +1,4 @@ -from pyprima import minimize +from pyprima import minimize, LinearConstraint as LC, NonlinearConstraint as NLC import numpy as np def test_eliminate_fixed_bounds(): @@ -13,3 +13,48 @@ def f(x): res = minimize(f, x0=np.array([1, 2, 3, 4, 5]), bounds=bounds) assert np.allclose(res.x, np.array([-0.5, -0.5, 1, 0, -0.5]), atol=1e-3) assert np.allclose(res.f, 1.75, atol=1e-3) + + +def test_eliminate_fixed_bounds_with_linear_constraints(): + # Ensure that the logic for fixed bounds also modifies linear constraints + # appropriately + + def f(x): + return np.sum(x**2) + + lb = [-1, None, None] + ub = [-1, None, None] + bounds = [(a, b) for a, b in zip(lb, ub)] + # Internally, the algorithm should modify lc to be a 1x2 matrix instead of 1x3, + # since it will modify x0 and the objective function to eliminate the first + # variable. If it fails to modify lc, we will get shape mismatch errors. + lc = LC(np.array([1, 1, 1]), lb=9, ub=15) + res = minimize(f, x0=np.array([1, 1, 1]), constraints=lc, bounds=bounds) + assert np.isclose(res.x[0], -1, atol=1e-6, rtol=1e-6) + assert np.isclose(res.x[1], 5, atol=1e-6, rtol=1e-6) + assert np.isclose(res.x[2], 5, atol=1e-6, rtol=1e-6) + assert np.isclose(res.f, 51, atol=1e-6, rtol=1e-6) + + +def test_eliminate_fixed_bounds_with_nonlinear_constraints(): + # Ensure that the logic for fixed bounds also modifies the nonlinear constraint + # function appropriately + + def f(x): + return np.sum(x**2) + + lb = [-1, None, None] + ub = [-1, None, None] + bounds = [(a, b) for a, b in zip(lb, ub)] + x0 = np.array([1, 1, 1]) + # Have the nonlinear constraint function operate on the last element of x, but be + # explicit about the length of x. This ensures that the test is still valid if the + # fixed bound is removed. If we simply used x[-1] this test would pass but it + # wouldn't actually test if we had properly modified the nonlinear constraint + # function after removing the fixed bounds + nlc = NLC(lambda x: x[len(x0)-1]**2, lb=9, ub=15) + res = minimize(f, x0=x0, constraints=nlc, bounds=bounds) + assert np.isclose(res.x[0], -1, atol=1e-6, rtol=1e-6) + assert np.isclose(res.x[1], 0, atol=1e-6, rtol=1e-6) + assert np.isclose(res.x[2], 3, atol=1e-6, rtol=1e-6) + assert np.isclose(res.f, 10, atol=1e-6, rtol=1e-6) \ No newline at end of file From e2e1fae51b54a12a3b14b5364702e18bed4fbb12 Mon Sep 17 00:00:00 2001 From: Nickolai Belakovski Date: Mon, 29 Sep 2025 11:09:38 -0400 Subject: [PATCH 2/5] Copy pyprima code for fixed constraints to python bindings Some minor modifications needed to be made. Notably pyprima's Result object uses 'f' for the function value whereas the bindings use 'fun'. Secondly, with the bindings we try to maintain the type that was provided, and so we made to make a minor modification as we changed the x0 that was provided to eliminate the fixed bounds. Ref: #250 --- python/prima/__init__.py | 63 ++++++++++++++++++++++++++++++++++++- python/prima/_common.py | 29 +++++++++++++++++ python/tests/test_bounds.py | 60 +++++++++++++++++++++++++++++++++++ 3 files changed, 151 insertions(+), 1 deletion(-) create mode 100644 python/tests/test_bounds.py diff --git a/python/prima/__init__.py b/python/prima/__init__.py index 33cbead755..0134755f02 100644 --- a/python/prima/__init__.py +++ b/python/prima/__init__.py @@ -10,6 +10,7 @@ from enum import Enum import numpy as np from ._common import _project +from ._common import get_arrays_tol class ConstraintType(Enum): @@ -115,6 +116,59 @@ def minimize(fun, x0, args=(), method=None, bounds=None, constraints=(), callbac 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], + ) + if isinstance(x0, np.ndarray): + x0 = np.array(x0)[~_fixed_idx] + elif x0_is_scalar: + raise Exception("You have provided a scalar x with fixed bounds, there is" \ + "no optimization to be done.") + else: + # In this case x is presumably a list, so we turn it into a numpy array + # for the convenience of indexing and then turn it back into a list + x0 = np.array(x0)[~_fixed_idx].tolist() + 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 + # Adjust linear_constraint + if linear_constraint: + new_lb = linear_constraint.lb - linear_constraint.A[:, _fixed_idx] @ _fixed_values + new_ub = linear_constraint.ub - linear_constraint.A[:, _fixed_idx] @ _fixed_values + new_A = linear_constraint.A[:, ~_fixed_idx] + linear_constraint = LinearConstraint(new_A, new_lb, new_ub) + # Adjust nonlinear constraints + if nonlinear_constraint_function: + original_nlc_fun = nonlinear_constraint_function + def fixed_nlc_fun(x): + newx = np.zeros(lenx0) + newx[_fixed_idx] = _fixed_values + newx[~_fixed_idx] = x + return original_nlc_fun(newx, *args) + nonlinear_constraint_function = fixed_nlc_fun + + # Project x0 onto the feasible set if nonlinear_constraint_function is None: result = _project(x0, lb, ub, {"linear": linear_constraint, "nonlinear": None}) @@ -142,7 +196,7 @@ def minimize(fun, x0, args=(), method=None, bounds=None, constraints=(), callbac options['nlconstr0'] = nlconstr0 options['m_nlcon'] = len(nlconstr0) - return _minimize( + result = _minimize( fun, x0, args, @@ -157,3 +211,10 @@ def minimize(fun, x0, args=(), method=None, bounds=None, constraints=(), callbac 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 diff --git a/python/prima/_common.py b/python/prima/_common.py index 950c247bf1..d94f1e8d0d 100644 --- a/python/prima/_common.py +++ b/python/prima/_common.py @@ -165,3 +165,32 @@ def _project(x0, lb, ub, constraints): return OptimizeResult(x=x0_c) return OptimizeResult(x=x0_c) + + +def get_arrays_tol(*arrays): + """ + Get a relative tolerance for a set of arrays. Borrowed from COBYQA + + Parameters + ---------- + *arrays: tuple + Set of `numpy.ndarray` to get the tolerance for. + + Returns + ------- + float + Relative tolerance for the set of arrays. + + Raises + ------ + ValueError + If no array is provided. + """ + if len(arrays) == 0: + raise ValueError("At least one array must be provided.") + size = max(array.size for array in arrays) + weight = max( + np.max(np.abs(array[np.isfinite(array)]), initial=1.0) + for array in arrays + ) + return 10.0 * eps * max(size, 1.0) * weight diff --git a/python/tests/test_bounds.py b/python/tests/test_bounds.py new file mode 100644 index 0000000000..9382376925 --- /dev/null +++ b/python/tests/test_bounds.py @@ -0,0 +1,60 @@ +from prima import minimize, LinearConstraint as LC, NonlinearConstraint as NLC +import numpy as np + +def test_eliminate_fixed_bounds(): + # Test the logic for detecting and eliminating fixed bounds + + def f(x): + return np.sum(x**2) + + lb = [-1, None, 1, None, -0.5] + ub = [-0.5, -0.5, None, None, -0.5] + bounds = [(a, b) for a, b in zip(lb, ub)] + res = minimize(f, x0=np.array([1, 2, 3, 4, 5]), bounds=bounds) + assert np.allclose(res.x, np.array([-0.5, -0.5, 1, 0, -0.5]), atol=1e-3) + assert np.allclose(res.fun, 1.75, atol=1e-3) + + +def test_eliminate_fixed_bounds_with_linear_constraints(): + # Ensure that the logic for fixed bounds also modifies linear constraints + # appropriately + + def f(x): + return np.sum(x**2) + + lb = [-1, None, None] + ub = [-1, None, None] + bounds = [(a, b) for a, b in zip(lb, ub)] + # Internally, the algorithm should modify lc to be a 1x2 matrix instead of 1x3, + # since it will modify x0 and the objective function to eliminate the first + # variable. If it fails to modify lc, we will get shape mismatch errors. + lc = LC(np.array([1, 1, 1]), lb=9, ub=15) + res = minimize(f, x0=np.array([1, 1, 1]), constraints=lc, bounds=bounds) + assert np.isclose(res.x[0], -1, atol=1e-6, rtol=1e-6) + assert np.isclose(res.x[1], 5, atol=1e-6, rtol=1e-6) + assert np.isclose(res.x[2], 5, atol=1e-6, rtol=1e-6) + assert np.isclose(res.fun, 51, atol=1e-6, rtol=1e-6) + + +def test_eliminate_fixed_bounds_with_nonlinear_constraints(): + # Ensure that the logic for fixed bounds also modifies the nonlinear constraint + # function appropriately + + def f(x): + return np.sum(x**2) + + lb = [-1, None, None] + ub = [-1, None, None] + bounds = [(a, b) for a, b in zip(lb, ub)] + x0 = np.array([1, 1, 1]) + # Have the nonlinear constraint function operate on the last element of x, but be + # explicit about the length of x. This ensures that the test is still valid if the + # fixed bound is removed. If we simply used x[-1] this test would pass but it + # wouldn't actually test if we had properly modified the nonlinear constraint + # function after removing the fixed bounds + nlc = NLC(lambda x: x[len(x0)-1]**2, lb=9, ub=15) + res = minimize(f, x0=x0, constraints=nlc, bounds=bounds) + assert np.isclose(res.x[0], -1, atol=1e-6, rtol=1e-6) + assert np.isclose(res.x[1], 0, atol=1e-6, rtol=1e-6) + assert np.isclose(res.x[2], 3, atol=1e-6, rtol=1e-6) + assert np.isclose(res.fun, 10, atol=1e-6, rtol=1e-6) \ No newline at end of file From 97433eec58bf5c40f102cc61509b3b6c11c8ac12 Mon Sep 17 00:00:00 2001 From: Nickolai Belakovski Date: Thu, 16 Oct 2025 16:55:29 -0400 Subject: [PATCH 3/5] Add code to deal with infeasible bounds in Python This was suggested in #251. --- .github/actions/spelling/allow.txt | 1 + pyprima/src/pyprima/common/_bounds.py | 5 +++++ pyprima/tests/test_bounds.py | 27 ++++++++++++++++++++++++++- python/prima/_bounds.py | 7 ++++++- python/tests/test_bounds.py | 27 ++++++++++++++++++++++++++- 5 files changed, 64 insertions(+), 3 deletions(-) diff --git a/.github/actions/spelling/allow.txt b/.github/actions/spelling/allow.txt index 96914224d6..f6294b5db3 100644 --- a/.github/actions/spelling/allow.txt +++ b/.github/actions/spelling/allow.txt @@ -2525,3 +2525,4 @@ Qoro HHmm MMdd ZAMB +excinfo diff --git a/pyprima/src/pyprima/common/_bounds.py b/pyprima/src/pyprima/common/_bounds.py index 89afe9529f..c5a81718ce 100644 --- a/pyprima/src/pyprima/common/_bounds.py +++ b/pyprima/src/pyprima/common/_bounds.py @@ -30,5 +30,10 @@ def process_bounds(bounds, lenx0): # If there were fewer bounds than variables, pad the rest with -/+ infinity lb = np.concatenate((lb, -np.inf*np.ones(lenx0 - len(lb)))) ub = np.concatenate((ub, np.inf*np.ones(lenx0 - len(ub)))) + + # Check the infeasibility of the bounds. + # TODO: Return a code instead of asserting because a feasibility check is a valid solution + infeasible = np.logical_not(lb <= ub) + assert not np.any(infeasible), f"Some of the provided bounds are infeasible. {infeasible=} {lb[infeasible]=}, {ub[infeasible]=}" return lb, ub diff --git a/pyprima/tests/test_bounds.py b/pyprima/tests/test_bounds.py index 80fda6a6d7..acccf235a4 100644 --- a/pyprima/tests/test_bounds.py +++ b/pyprima/tests/test_bounds.py @@ -1,5 +1,6 @@ from pyprima import minimize, LinearConstraint as LC, NonlinearConstraint as NLC import numpy as np +import pytest def test_eliminate_fixed_bounds(): # Test the logic for detecting and eliminating fixed bounds @@ -57,4 +58,28 @@ def f(x): assert np.isclose(res.x[0], -1, atol=1e-6, rtol=1e-6) assert np.isclose(res.x[1], 0, atol=1e-6, rtol=1e-6) assert np.isclose(res.x[2], 3, atol=1e-6, rtol=1e-6) - assert np.isclose(res.f, 10, atol=1e-6, rtol=1e-6) \ No newline at end of file + assert np.isclose(res.f, 10, atol=1e-6, rtol=1e-6) + + +def test_infeasible_bounds(): + def f(x): + return np.sum(x**2) + + lb = [1, None, None] + ub = [-1, None, None] + bounds = [(a, b) for a, b in zip(lb, ub)] + with pytest.raises(AssertionError) as excinfo: + minimize(f, x0=np.array([1, 2, 3]), bounds=bounds) + assert str(excinfo.value) == "Some of the provided bounds are infeasible. infeasible=array([ True, False, False]) lb[infeasible]=array([1.]), ub[infeasible]=array([-1.])" + + +def test_infeasible_bounds_nan(): + def f(x): + return np.sum(x**2) + + lb = [np.nan, None, None] + ub = [-1, None, None] + bounds = [(a, b) for a, b in zip(lb, ub)] + with pytest.raises(AssertionError) as excinfo: + minimize(f, x0=np.array([1, 2, 3]), bounds=bounds) + assert str(excinfo.value) == "Some of the provided bounds are infeasible. infeasible=array([ True, False, False]) lb[infeasible]=array([nan]), ub[infeasible]=array([-1.])" diff --git a/python/prima/_bounds.py b/python/prima/_bounds.py index 89afe9529f..a917a0bc1e 100644 --- a/python/prima/_bounds.py +++ b/python/prima/_bounds.py @@ -30,5 +30,10 @@ def process_bounds(bounds, lenx0): # If there were fewer bounds than variables, pad the rest with -/+ infinity lb = np.concatenate((lb, -np.inf*np.ones(lenx0 - len(lb)))) ub = np.concatenate((ub, np.inf*np.ones(lenx0 - len(ub)))) - + + # Check the infeasibility of the bounds. + # TODO: Return a code instead of asserting because a feasibility check is a valid solution + infeasible = np.logical_not(lb <= ub) + assert not np.any(infeasible), f"Some of the provided bounds are infeasible. {infeasible=} {lb[infeasible]=}, {ub[infeasible]=}" + return lb, ub diff --git a/python/tests/test_bounds.py b/python/tests/test_bounds.py index 9382376925..f4153ef9ff 100644 --- a/python/tests/test_bounds.py +++ b/python/tests/test_bounds.py @@ -1,5 +1,6 @@ from prima import minimize, LinearConstraint as LC, NonlinearConstraint as NLC import numpy as np +import pytest def test_eliminate_fixed_bounds(): # Test the logic for detecting and eliminating fixed bounds @@ -57,4 +58,28 @@ def f(x): assert np.isclose(res.x[0], -1, atol=1e-6, rtol=1e-6) assert np.isclose(res.x[1], 0, atol=1e-6, rtol=1e-6) assert np.isclose(res.x[2], 3, atol=1e-6, rtol=1e-6) - assert np.isclose(res.fun, 10, atol=1e-6, rtol=1e-6) \ No newline at end of file + assert np.isclose(res.fun, 10, atol=1e-6, rtol=1e-6) + + +def test_infeasible_bounds(): + def f(x): + return np.sum(x**2) + + lb = [1, None, None] + ub = [-1, None, None] + bounds = [(a, b) for a, b in zip(lb, ub)] + with pytest.raises(AssertionError) as excinfo: + minimize(f, x0=np.array([1, 2, 3]), bounds=bounds) + assert str(excinfo.value) == "Some of the provided bounds are infeasible. infeasible=array([ True, False, False]) lb[infeasible]=array([1.]), ub[infeasible]=array([-1.])" + + +def test_infeasible_bounds_nan(): + def f(x): + return np.sum(x**2) + + lb = [np.nan, None, None] + ub = [-1, None, None] + bounds = [(a, b) for a, b in zip(lb, ub)] + with pytest.raises(AssertionError) as excinfo: + minimize(f, x0=np.array([1, 2, 3]), bounds=bounds) + assert str(excinfo.value) == "Some of the provided bounds are infeasible. infeasible=array([ True, False, False]) lb[infeasible]=array([nan]), ub[infeasible]=array([-1.])" From affd6d34fca58fa2f0a21f33f97462d483a238dc Mon Sep 17 00:00:00 2001 From: Nickolai Belakovski Date: Tue, 6 Jan 2026 18:24:55 +0800 Subject: [PATCH 4/5] Return early in the case of all fixed x Also evalaute the objective and nonlinear constraint functions A 'bug' in pyprima was found and fixed, in which it was returning both linear and nonlinear constraint evaluations. It has been fixed to return only nonlinear constraint evaluations, to be consistent with the Fortran implementation. --- pyprima/src/pyprima/__init__.py | 30 ++++++++++++++-- pyprima/src/pyprima/cobyla/cobyla.py | 2 ++ pyprima/src/pyprima/common/infos.py | 1 + pyprima/tests/test_bounds.py | 51 ++++++++++++++++++++++++++++ python/_prima.cpp | 3 ++ python/prima/__init__.py | 37 ++++++++++++++++++-- python/prima/infos.py | 5 +++ python/tests/test_bounds.py | 51 ++++++++++++++++++++++++++++ 8 files changed, 175 insertions(+), 5 deletions(-) create mode 100644 python/prima/infos.py diff --git a/pyprima/src/pyprima/__init__.py b/pyprima/src/pyprima/__init__.py index 1c7f54160c..2feb40ac80 100644 --- a/pyprima/src/pyprima/__init__.py +++ b/pyprima/src/pyprima/__init__.py @@ -8,8 +8,9 @@ from .common._bounds import process_bounds from enum import Enum from .common._project import _project +from .common.infos import FIXED_SUCCESS from .common.linalg import get_arrays_tol -from .cobyla.cobyla import cobyla +from .cobyla.cobyla import cobyla, COBYLAResult import numpy as np from collections.abc import Iterable @@ -132,7 +133,7 @@ def scalar_fun(x): (lb <= ub) & (np.abs(lb - ub) < tol) ) - if any(_fixed_idx): + if any(_fixed_idx) and not all(_fixed_idx): _fixed_values = 0.5 * ( lb[_fixed_idx] + ub[_fixed_idx] ) @@ -194,6 +195,31 @@ def calcfc(x): constr = np.zeros(0) return f, constr + if all(_fixed_idx): + x=0.5 * ( lb[_fixed_idx] + ub[_fixed_idx] ) + x = np.clip( + x, + lb[_fixed_idx], + ub[_fixed_idx], + ) + f, nlconstr = calcfc(x) + cstrv = max(max((A_ineq @ x) - b_ineq) if A_ineq is not None else 0, + max((abs((A_eq @ x) - b_eq))) if A_eq is not None else 0, + max(np.append(0, nlconstr))) + result = COBYLAResult( + x=x, + f=f, + constr=nlconstr, + cstrv=cstrv, + nf=1, + xhist=[x], + fhist=[f], + chist=[cstrv], + conhist=[nlconstr], + info=FIXED_SUCCESS + ) + return result + f0, nlconstr0 = calcfc(x0) if 'quiet' in options: diff --git a/pyprima/src/pyprima/cobyla/cobyla.py b/pyprima/src/pyprima/cobyla/cobyla.py index 11772a249e..e9f190668b 100644 --- a/pyprima/src/pyprima/cobyla/cobyla.py +++ b/pyprima/src/pyprima/cobyla/cobyla.py @@ -481,6 +481,8 @@ def calcfc(x): callback ) + constr = constr[mmm - m_nlcon:] + return COBYLAResult(x, f, constr, cstrv, nf, xhist, fhist, chist, conhist, info) diff --git a/pyprima/src/pyprima/common/infos.py b/pyprima/src/pyprima/common/infos.py index 27d0f1e144..dc2975576b 100644 --- a/pyprima/src/pyprima/common/infos.py +++ b/pyprima/src/pyprima/common/infos.py @@ -13,6 +13,7 @@ FTARGET_ACHIEVED = 1 TRSUBP_FAILED = 2 MAXFUN_REACHED = 3 +FIXED_SUCCESS = 13 MAXTR_REACHED = 20 NAN_INF_X = -1 NAN_INF_F = -2 diff --git a/pyprima/tests/test_bounds.py b/pyprima/tests/test_bounds.py index acccf235a4..1c4f8b3a25 100644 --- a/pyprima/tests/test_bounds.py +++ b/pyprima/tests/test_bounds.py @@ -83,3 +83,54 @@ def f(x): with pytest.raises(AssertionError) as excinfo: minimize(f, x0=np.array([1, 2, 3]), bounds=bounds) assert str(excinfo.value) == "Some of the provided bounds are infeasible. infeasible=array([ True, False, False]) lb[infeasible]=array([nan]), ub[infeasible]=array([-1.])" + + +@pytest.mark.parametrize('with_Aineq', [False, True]) +@pytest.mark.parametrize('with_Aeq', [False, True]) +@pytest.mark.parametrize('with_nlcon', [False, True]) +def test_all_fixed(with_nlcon, with_Aeq, with_Aineq): + def f(x): + return np.sum(x**2) + + lb = [-1, 2, 4] + ub = [-1, 2, 4] + bounds = [(a, b) for a, b in zip(lb, ub)] + nlc = NLC(lambda x: x[2]**2, lb=-np.inf, ub=0) + lc_eq = LC([0, 0, 1], lb=21, ub=21) + lc_ineq = LC([0, 0, 1], lb=22, ub=23) + constraints = [] + if with_nlcon: + constraints.append(nlc) + if with_Aeq: + constraints.append(lc_eq) + if with_Aineq: + constraints.append(lc_ineq) + + result = minimize(f, x0=np.array([1, 2, 3]), bounds=bounds, constraints=constraints) + assert all(result.x == [-1, 2, 4]) + assert result.f == 21 + assert result.nf == 1 + if not with_nlcon and not with_Aeq and not with_Aineq: + assert np.array_equal(result.constr, []) + assert result.cstrv == 0 + if not with_nlcon and not with_Aeq and with_Aineq: + assert np.array_equal(result.constr, []) + assert result.cstrv == 18 + if not with_nlcon and with_Aeq and not with_Aineq: + assert np.array_equal(result.constr, []) + assert result.cstrv == 17 + if not with_nlcon and with_Aeq and with_Aineq: + assert np.array_equal(result.constr, []) + assert result.cstrv == 18 + if with_nlcon and not with_Aeq and not with_Aineq: + assert np.array_equal(result.constr, [16]) + assert result.cstrv == 16 + if with_nlcon and not with_Aeq and with_Aineq: + assert np.array_equal(result.constr, [16]) + assert result.cstrv == 18 + if with_nlcon and with_Aeq and not with_Aineq: + assert np.array_equal(result.constr, [16]) + assert result.cstrv == 17 + if with_nlcon and with_Aeq and with_Aineq: + assert np.array_equal(result.constr, [16]) + assert result.cstrv == 18 diff --git a/python/_prima.cpp b/python/_prima.cpp index 3dda286693..349c794c2a 100644 --- a/python/_prima.cpp +++ b/python/_prima.cpp @@ -29,6 +29,8 @@ class SelfCleaningPyObject { }; struct PRIMAResult { + // Blank constructor for construction from Python + PRIMAResult() {} // Construct PRIMAResult from prima_result_t PRIMAResult(const prima_result_t& result, const int num_vars, const int num_constraints, const std::string method) : x(num_vars, result.x), @@ -89,6 +91,7 @@ PYBIND11_MODULE(_prima, m) { #endif py::class_(m, "PRIMAResult") + .def(py::init<>()) .def_readwrite("x", &PRIMAResult::x) .def_readwrite("success", &PRIMAResult::success) .def_readwrite("status", &PRIMAResult::status) diff --git a/python/prima/__init__.py b/python/prima/__init__.py index 0134755f02..119df48db1 100644 --- a/python/prima/__init__.py +++ b/python/prima/__init__.py @@ -1,4 +1,4 @@ -from ._prima import minimize as _minimize, __version__, PRIMAMessage +from ._prima import minimize as _minimize, __version__, PRIMAMessage, PRIMAResult # 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 ._nonlinear_constraints import process_nl_constraints @@ -11,6 +11,7 @@ import numpy as np from ._common import _project from ._common import get_arrays_tol +from .infos import FIXED_SUCCESS class ConstraintType(Enum): @@ -125,7 +126,7 @@ def minimize(fun, x0, args=(), method=None, bounds=None, constraints=(), callbac (lb <= ub) & (np.abs(lb - ub) < tol) ) - if any(_fixed_idx): + if any(_fixed_idx) and not all(_fixed_idx): _fixed_values = 0.5 * ( lb[_fixed_idx] + ub[_fixed_idx] ) @@ -183,7 +184,7 @@ def fixed_nlc_fun(x): else: A_eq, b_eq, A_ineq, b_ineq = None, None, None, None - if nonlinear_constraint_function is not None: + if nonlinear_constraint_function is not None and not all(_fixed_idx): # 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 @@ -196,6 +197,36 @@ def fixed_nlc_fun(x): options['nlconstr0'] = nlconstr0 options['m_nlcon'] = len(nlconstr0) + if all(_fixed_idx): + x = 0.5 * ( + lb[_fixed_idx] + ub[_fixed_idx] + ) + x = np.clip( + x, + lb[_fixed_idx], + ub[_fixed_idx], + ) + success = True + status = FIXED_SUCCESS + message = "All variables were fixed by the provided bounds." + fun = fun(x) + nfev = 1 + nlconstr = nonlinear_constraint_function(x) if nonlinear_constraint_function is not None else np.array([]) + maxcv = max(max((A_ineq @ x) - b_ineq) if A_ineq is not None else 0, + max((abs((A_eq @ x) - b_eq))) if A_eq is not None else 0, + max(np.append(0, nlconstr))) + result = PRIMAResult() + result.x = x + result.success = success + result.status = status + result.message = message + result.fun = fun + result.nfev = nfev + result.maxcv = maxcv + result.nlconstr = nlconstr + result.method = method + return result + result = _minimize( fun, x0, diff --git a/python/prima/infos.py b/python/prima/infos.py new file mode 100644 index 0000000000..77ffea2ca1 --- /dev/null +++ b/python/prima/infos.py @@ -0,0 +1,5 @@ +''' +This module will be merged with pyprima infos.py in the near future +''' + +FIXED_SUCCESS = 13 \ No newline at end of file diff --git a/python/tests/test_bounds.py b/python/tests/test_bounds.py index f4153ef9ff..241223d312 100644 --- a/python/tests/test_bounds.py +++ b/python/tests/test_bounds.py @@ -83,3 +83,54 @@ def f(x): with pytest.raises(AssertionError) as excinfo: minimize(f, x0=np.array([1, 2, 3]), bounds=bounds) assert str(excinfo.value) == "Some of the provided bounds are infeasible. infeasible=array([ True, False, False]) lb[infeasible]=array([nan]), ub[infeasible]=array([-1.])" + + +@pytest.mark.parametrize('with_Aineq', [False, True]) +@pytest.mark.parametrize('with_Aeq', [False, True]) +@pytest.mark.parametrize('with_nlcon', [False, True]) +def test_all_fixed(with_nlcon, with_Aeq, with_Aineq): + def f(x): + return np.sum(x**2) + + lb = [-1, 2, 4] + ub = [-1, 2, 4] + bounds = [(a, b) for a, b in zip(lb, ub)] + nlc = NLC(lambda x: x[2]**2, lb=-np.inf, ub=0) + lc_eq = LC([0, 0, 1], lb=21, ub=21) + lc_ineq = LC([0, 0, 1], lb=22, ub=23) + constraints = [] + if with_nlcon: + constraints.append(nlc) + if with_Aeq: + constraints.append(lc_eq) + if with_Aineq: + constraints.append(lc_ineq) + + result = minimize(f, x0=np.array([1, 2, 3]), bounds=bounds, constraints=constraints) + assert all(result.x == [-1, 2, 4]) + assert result.fun == 21 + assert result.nfev == 1 + if not with_nlcon and not with_Aeq and not with_Aineq: + assert np.array_equal(result.nlconstr, []) + assert result.maxcv == 0 + if not with_nlcon and not with_Aeq and with_Aineq: + assert np.array_equal(result.nlconstr, []) + assert result.maxcv == 18 + if not with_nlcon and with_Aeq and not with_Aineq: + assert np.array_equal(result.nlconstr, []) + assert result.maxcv == 17 + if not with_nlcon and with_Aeq and with_Aineq: + assert np.array_equal(result.nlconstr, []) + assert result.maxcv == 18 + if with_nlcon and not with_Aeq and not with_Aineq: + assert np.array_equal(result.nlconstr, [16]) + assert result.maxcv == 16 + if with_nlcon and not with_Aeq and with_Aineq: + assert np.array_equal(result.nlconstr, [16]) + assert result.maxcv == 18 + if with_nlcon and with_Aeq and not with_Aineq: + assert np.array_equal(result.nlconstr, [16]) + assert result.maxcv == 17 + if with_nlcon and with_Aeq and with_Aineq: + assert np.array_equal(result.nlconstr, [16]) + assert result.maxcv == 18 From f5e17e5c4a745fcbb3e380b8b8a88a08a6116880 Mon Sep 17 00:00:00 2001 From: Nickolai Belakovski Date: Wed, 7 Jan 2026 16:18:48 +0800 Subject: [PATCH 5/5] Minor PR feedback --- pyprima/src/pyprima/__init__.py | 7 +++++-- python/prima/__init__.py | 10 +++++----- 2 files changed, 10 insertions(+), 7 deletions(-) diff --git a/pyprima/src/pyprima/__init__.py b/pyprima/src/pyprima/__init__.py index 2feb40ac80..704fbff75a 100644 --- a/pyprima/src/pyprima/__init__.py +++ b/pyprima/src/pyprima/__init__.py @@ -131,9 +131,12 @@ def scalar_fun(x): tol = get_arrays_tol(lb, ub) _fixed_idx = ( (lb <= ub) - & (np.abs(lb - ub) < tol) + & (ub <= lb + tol) ) if any(_fixed_idx) and not all(_fixed_idx): + # We should NOT reduce the problem if all variables are fixed. Otherwise, Aineq would be [], and + # then bineq will be set to [] in the end. In this way, we lose completely the information in + # these constraints. Consequently, we cannot evaluate the constraint violation correctly when needed. _fixed_values = 0.5 * ( lb[_fixed_idx] + ub[_fixed_idx] ) @@ -246,7 +249,7 @@ def calcfc(x): ) if any(_fixed_idx): - newx = np.zeros(lenx0) + newx = np.zeros(lenx0) + np.nan newx[_fixed_idx] = _fixed_values newx[~_fixed_idx] = result.x result.x = newx diff --git a/python/prima/__init__.py b/python/prima/__init__.py index 119df48db1..a6eaf6ec78 100644 --- a/python/prima/__init__.py +++ b/python/prima/__init__.py @@ -124,9 +124,12 @@ def minimize(fun, x0, args=(), method=None, bounds=None, constraints=(), callbac tol = get_arrays_tol(lb, ub) _fixed_idx = ( (lb <= ub) - & (np.abs(lb - ub) < tol) + & (ub <= lb + tol) ) if any(_fixed_idx) and not all(_fixed_idx): + # We should NOT reduce the problem if all variables are fixed. Otherwise, Aineq would be [], and + # then bineq will be set to [] in the end. In this way, we lose completely the information in + # these constraints. Consequently, we cannot evaluate the constraint violation correctly when needed. _fixed_values = 0.5 * ( lb[_fixed_idx] + ub[_fixed_idx] ) @@ -137,9 +140,6 @@ def minimize(fun, x0, args=(), method=None, bounds=None, constraints=(), callbac ) if isinstance(x0, np.ndarray): x0 = np.array(x0)[~_fixed_idx] - elif x0_is_scalar: - raise Exception("You have provided a scalar x with fixed bounds, there is" \ - "no optimization to be done.") else: # In this case x is presumably a list, so we turn it into a numpy array # for the convenience of indexing and then turn it back into a list @@ -244,7 +244,7 @@ def fixed_nlc_fun(x): ) if any(_fixed_idx): - newx = np.zeros(lenx0) + newx = np.zeros(lenx0) + np.nan newx[_fixed_idx] = _fixed_values newx[~_fixed_idx] = result.x result.x = newx