Make pyprima update constraints when eliminating fixed variables and bring the code to python bindings

.github/actions/spelling/allow.txt

@@ -2525,3 +2525,4 @@ Qoro
 HHmm
 MMdd
 ZAMB
+excinfo

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
 
@@ -130,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):
+    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]
         )
@@ -151,6 +155,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
@@ -179,6 +198,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:
@@ -205,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

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)
 
 

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

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

pyprima/tests/test_bounds.py

@@ -1,5 +1,6 @@
-from pyprima import minimize
+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
@@ -13,3 +14,123 @@ 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)
+
+
+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.])"
+
+
+@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

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_<PRIMAResult>(m, "PRIMAResult")
+      .def(py::init<>())
       .def_readwrite("x", &PRIMAResult::x)
       .def_readwrite("success", &PRIMAResult::success)
       .def_readwrite("status", &PRIMAResult::status)