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

Author: nbelakovski

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

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

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)