Refactor the processing of nonlinear constraints

.github/actions/spelling/allow.txt

@@ -2248,3 +2248,4 @@ pycutest
 excludelist
 slsqp
 optiprofiler
+manylinux

python/prima/__init__.py

@@ -32,10 +32,10 @@ def get_constraint_type(constraint):
         raise ValueError("Constraint type not recognized")
 
 
-def process_constraints(constraints):
+def process_constraints(constraints, x0):
     # First throw it back if it's an empty tuple
     if not constraints:
-        return None, None
+        return None, 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
@@ -59,9 +59,10 @@ def process_constraints(constraints):
             raise ValueError("Constraint type not recognized")
 
     if len(nonlinear_constraints) > 0:
-        nonlinear_constraint_function = process_nl_constraints(nonlinear_constraints)
+        nonlinear_constraint_function, nlconstr0 = process_nl_constraints(nonlinear_constraints, x0)
     else:
         nonlinear_constraint_function = None
+        nlconstr0 = None
 
     # Determine if we have multiple linear constraints, just 1, or none, and process accordingly
     if len(linear_constraints) > 1:
@@ -71,12 +72,12 @@ def process_constraints(constraints):
     else:
         linear_constraint = None
 
-    return linear_constraint, nonlinear_constraint_function
+    return linear_constraint, nonlinear_constraint_function, nlconstr0
 
 
 def minimize(fun, x0, args=(), method=None, bounds=None, constraints=(), callback=None, options=None):
 
-    linear_constraint, nonlinear_constraint_function = process_constraints(constraints)
+    linear_constraint, nonlinear_constraint_function, nlconstr0 = process_constraints(constraints, x0)
 
     quiet = options.get("quiet", True) if options is not None else True
 
@@ -129,11 +130,9 @@ def minimize(fun, x0, args=(), method=None, bounds=None, constraints=(), callbac
         A_eq, b_eq, A_ineq, b_ineq = None, None, None, None
 
     if nonlinear_constraint_function is not None:
-        nlconstr0 = nonlinear_constraint_function(x0)
-        m_nlcon = len(nlconstr0) if hasattr(nlconstr0, "__len__") else 1
+        m_nlcon = len(nlconstr0)
         f0 = fun(x0)
     else:
-        nlconstr0 = None
         m_nlcon = None
         f0 = None
 

python/prima/_nonlinear_constraints.py

@@ -1,100 +1,79 @@
 import numpy as np
 
-def process_single_nl_constraint(nlc):
+def process_single_nl_constraint(nlc, x0):
     '''
     The Python interfaces receives the constraints as lb <= constraint(x) <= ub, 
-    but the Fortran backend expects the nonlinear constraints to be constr(x) <= 0.
+    but the Fortran backend expects the nonlinear constraints to be constraint(x) <= 0.
     Thus a conversion is needed.
     
-    This function will take a NonlinearConstraint and return a new function which incorporates
-    the lower and upper bounds in such a way as to satisfy the expectation of the Fortran
-    backend. This new function first calls the original function from the provided NonlinearConstraint
-    and obtains the values, and then passes those values to a second function which combines
-    them with the lower and upper bounds as appropriate. The main work here is creating that
-    second function.
+    In addition to the conversion, we would like to check the size of the output of the
+    constraint function as compared to the provided lb and ub, and so we must run the
+    constraint function here. Since we presume it is expensive to run the constraint function,
+    we must also return the output of the constraint function to the caller, so that they may
+    avoid running the constraint function unnecessarily.
     
-    In general to convert lb <= constraint(x) <= ub to newconstraint(x) <= 0, all we need to do
-    is define newconstraint(x) as [constraint(x) - ub, lb - constraint(x)]. There are some details
-    regarding the type and content of lb and ub as detailed below.
+    In order to accomplish all these goals, we first run the constraint function, then we
+    upgrade lb/ub to vectors if they are not already (while checking their sizes in the process
+    and raising exceptions if necessary), and then we create a function to transform the output
+    of the constraint function to the form expected by the Fortran backend.
+
+    '''
 
+    # Run the constraint function and upgrade the lb/ub to the correct size,
+    # while also checking that if they were already vectors that they had the correct
+    # size in the first place.
+    nlconstr0 = nlc.fun(x0)
+    nlconstr0 = np.atleast_1d(np.array(nlconstr0, dtype=np.float64))
+    lb, ub = _upgrade_lb_ub_to_vectors(nlc.lb, nlc.ub, nlconstr0)
 
-    The upper and lower bounds could be either scalars or vectors, and so we have 4 possible
-    cases to take into account:
+    # Check if any bounds contain -inf and +inf simultaneously
+    if np.any((lb == -np.inf) & (ub == np.inf)):
+        raise ValueError("A NonlinearConstraint was provided without specifying lower or upper bounds")
 
-    1. Both lb and ub are scalars
-        In this case we have a further 4 cases!
-        
-        1a. lb == -inf and ub == inf
-            This case makes no sense since it means that there are effectively no constraints,
-            so we raise an exception.
-        1b. lb == -inf and ub != inf
-            This is a one sided constraint, so we can define newconstraint(x) as constraint(x) - ub
-        1c. lb != -inf and ub == inf
-            This is a one sided constraint, so we can define newconstraint(x) as lb - constraint(x)
-        1d. lb != -inf and ub != inf
-            Since we have both constraints we define newconstraint(x) as [constraint(x) - ub, lb - constraint(x)]
-    2. lb is a scalar and ub is a vector
-        In this case we have a further 2 cases.
-        
-        2a. lb == -inf
-            This is a one sided constraint, so we can define newconstraint(x) as [constraint(x) - ub], however we
-            should first check if there is any inf in ub and raise an exception since this would be same case as 1a.
-        2b. lb != -inf
-            In this case we can have inf in ub, so newconstraint needs to check for that and omit it from the constraints.
-            See the code for the exact method of accomplishing this.
-    3. lb is a vector and ub is a scalar
-        This is the same as case 2, but with lb and ub reversed.
-    4. Both lb and ub are vectors
-        There are no subcases here, other than checking for -np.inf in lb and np.inf in ub and removing those constraints.
-        However we also check for any particular index where ub is inf and lb is -inf simultaneously and raise an exception
-        if so, since again this is like case 1a.
-
-    '''
-    lb_is_scalar = not hasattr(nlc.lb, "__len__")
-    ub_is_scalar = not hasattr(nlc.ub, "__len__")
-    if lb_is_scalar and ub_is_scalar:
-        if nlc.lb == -np.inf and nlc.ub == np.inf:
-            raise ValueError("A NonlinearConstraint was provided without specifying lower or upper bounds")
-        elif nlc.lb == -np.inf:
-            align_constraint_values = lambda values: values - nlc.ub
-        elif nlc.ub == np.inf:
-            align_constraint_values = lambda values: nlc.lb - values
-        else:
-            align_constraint_values = lambda values: np.concatenate((values - nlc.ub, nlc.lb - values))
-    elif lb_is_scalar and not ub_is_scalar:
-        if nlc.lb == -np.inf:
-            if np.inf in nlc.ub:
-                raise ValueError("A NonlinearConstraint was provided without specifying lower or upper bounds")
-            align_constraint_values = lambda values: values - nlc.ub
-        else:
-            align_constraint_values = lambda values: np.concatenate(([vi - ub_ii for ub_ii, vi in zip(nlc.ub, values) if ub_ii < np.inf],
-                                        nlc.lb - values))
-    elif not lb_is_scalar and ub_is_scalar:
-        if nlc.ub == np.inf:
-            if -np.inf in nlc.lb:
-                raise ValueError("A NonlinearConstraint was provided without specifying lower or upper bounds")
-            align_constraint_values = lambda values: nlc.lb - values
-        else:
-            align_constraint_values = lambda values: np.concatenate((values - nlc.ub,
-                    [lb_ii - vi for lb_ii, vi in zip(nlc.lb, values) if lb_ii > -np.inf]))
-    else:
-        if np.any((nlc.lb == -np.inf) & (nlc.ub == np.inf)):
-            raise ValueError("A NonlinearConstraint was provided without specifying lower or upper bounds")
-        align_constraint_values = lambda values: np.concatenate(([vi - ub_ii for ub_ii, vi in zip(nlc.ub, values) if ub_ii < np.inf],
-                                   [lb_ii - vi for lb_ii, vi in zip(nlc.lb, values) if lb_ii > -np.inf]))
+    align_constraint_values = lambda values: np.concatenate(([vi - ub_ii for ub_ii, vi in zip(ub, values) if ub_ii < np.inf],
+                                [lb_ii - vi for lb_ii, vi in zip(lb, values) if lb_ii > -np.inf]))
     def newconstraint(x):
         values = np.atleast_1d(np.array(nlc.fun(x), dtype=np.float64))
         return align_constraint_values(values)
-    return newconstraint
+    nlconstr0 = align_constraint_values(nlconstr0)
+    return newconstraint, nlconstr0
 
 
-def process_nl_constraints(nlcs):
+def process_nl_constraints(nlcs, x0):
     functions = []
+    nlconstr0 = np.empty(0)
     for nlc in nlcs:
-        functions.append(process_single_nl_constraint(nlc))
+        fun_i, nlconstr0_i = process_single_nl_constraint(nlc, x0)
+        functions.append(fun_i)
+        nlconstr0 = np.concatenate((nlconstr0, nlconstr0_i))
     def constraint_function(x):
         values = np.empty(0)
         for fun in functions:
             values = np.concatenate((values, fun(x)))
         return values
-    return constraint_function
+    return constraint_function, nlconstr0
+
+
+def _upgrade_lb_ub_to_vectors(lb, ub, nlconstr0):
+    '''
+    Make sure length of lb/ub align with length of nlconstr0
+    '''
+    lb_is_scalar = not hasattr(lb, "__len__")
+    ub_is_scalar = not hasattr(ub, "__len__")
+    len_nlconstr0 = len(nlconstr0)
+    if lb_is_scalar and ub_is_scalar:
+        return np.array([lb]*len_nlconstr0, dtype=np.float64), np.array([ub]*len_nlconstr0, dtype=np.float64)
+    elif lb_is_scalar and not ub_is_scalar:
+        if len(ub) != len_nlconstr0:
+            raise ValueError("The number of elements in the constraint function's output does not match the number of elements in the upper bound.")
+        return np.array([lb]*len_nlconstr0, dtype=np.float64), np.array(ub, dtype=np.float64)
+    elif not lb_is_scalar and ub_is_scalar:
+        if len(lb) != len_nlconstr0:
+            raise ValueError("The number of elements in the constraint function's output does not match the number of elements in the lower bound.")
+        return np.array(lb, dtype=np.float64), np.array([ub]*len_nlconstr0, dtype=np.float64)
+    else:
+        if len(lb) != len_nlconstr0:
+            raise ValueError("The number of elements in the constraint function's output does not match the number of elements in the lower bound.")
+        if len(ub) != len_nlconstr0:
+            raise ValueError("The number of elements in the constraint function's output does not match the number of elements in the upper bound.")
+        return np.array(lb, dtype=np.float64), np.array(ub, dtype=np.float64)

python/tests/README.md

@@ -61,8 +61,18 @@ order to run that test by itself.
 - [x] compatible with PDFO API - test_compatibility_pdfo.py::test_pdfo
 
 
-## Requirements for processing of lists of constraints
-- [x] providing a list of nonlinear constraint functions without providing either the total dimension or the dimension of each function successfully determines the total number of constraints - test_process_nonlinear_constraints.py::test_multiple_nl_constraints_various_data_types
+## Requirements for processing nonlinear constraints
+- [x] providing a list of nonlinear constraints with either scalar or vector for lb/ub correctly constructs the constraint function - test_process_nonlinear_constraints.py::test_multiple_nl_constraints_various_data_types
+- [x] providing a single nonlinear constraints with either scalar or vector for lb/ub correctly constructs the constraint function - test_process_nonlinear_constraints.py::test_single_nl_constraint
+- [x] providing a nonlinear constraint with upper and/or lower bounds as vectors of different length than the length of the output of the constraint function raises an exception - test_process_nonlinear_constraints.py::test_length_nlc_fun_not_equal_to_length_lb_ub
+- [x] providing a nonlinear constraint with upper and lower bound as vectors of where upper bounds has correct length and lower bounds does not as compared to the length of the output of the constraint function raises an exception - test_process_nonlinear_constraints.py::test_length_nlc_fun_ne_to_length_ub
+- [x] providing a nonlinear constraint with scalar lb at -inf and scalar ub at inf raises an exception - test_process_nonlinear_constraints.py::test_lb_neg_inf_ub_inf_raises
+- [x] providing a nonlinear constraint with scalar lb at -inf and vector ub containing at least 1 inf raises an exception - test_process_nonlinear_constraints.py::test_lb_neg_inf_ub_vector_w_inf_raises
+- [x] providing a nonlinear constraint with vector lb containing at least 1 -inf and scalar ub at inf raises an exception - test_process_nonlinear_constraints.py::test_lb_vector_with_neg_inf_ub_inf_raises
+- [x] providing a nonlinear constraint with vector lb and vector ub containing -inf and +inf at the same index raises an exception - test_process_nonlinear_constraints.py::test_lb_vector_with_neg_inf_ub_vector_w_inf_at_same_index_raises
+
+
+## Requirements for processing linear constraints
 - [x] providing a list of linear constraints leads to them being successfully combined - test_process_linear_constraints.py::test_multiple_linear_constraints_implementation
 - [x] providing a list of linear constraints leads to them being successfully combined and applied - test_process_linear_constraints.py::test_multiple_linear_constraints_high_level
 - [x] providing linear constraints that contain both equality and inequality constraints leads to a correct separation into A_eq/b_eq and A_ineq/b_ineq - test_process_linear_constraints.py::test_separate_LC_into_eq_and_ineq
@@ -78,7 +88,7 @@ order to run that test by itself.
 - [x] providing an unsupported constraint type raises an exception - test_constraint_types.py::test_unsupported_type_raises_exception
 
 
-## Requirements for warnings related to method selection
+## Requirements for warnings and exceptions related to method selection
 - [x] providing a method that is not one of the 5 provided raises an exception - test_method_selection_warnings_and_exceptions.py::test_method_not_recognized
 - [x] providing nonlinear constraints to a method that cannot handle them raises an exception - test_method_selection_warnings_and_exceptions.py::test_providing_nonlinear_constraints_to_non_cobyla_method
 - [x] providing linear constraints to a method that cannot handle them raises an exception - test_method_selection_warnings_and_exceptions.py::test_providing_linear_constraints_to_non_cobyla_non_lincoa_method

python/tests/test_process_nonlinear_constraints.py

@@ -1,5 +1,5 @@
 import numpy as np
-from prima import NonlinearConstraint, process_nl_constraints
+from prima import NonlinearConstraint, process_nl_constraints, minimize
 import pytest
 
 
@@ -11,14 +11,68 @@ def test_multiple_nl_constraints_various_data_types(lb1, ub1, lb2, ub2):
     nlc1 = NonlinearConstraint(lambda x: x, lb=lb1, ub=ub1)
     nlc2 = NonlinearConstraint(lambda x: [x, x], lb=lb2, ub=ub2)
     nlcs = [nlc1, nlc2]
-    constr_func = process_nl_constraints(nlcs)
+    x0 = 1
+    constr_func, nlconstr0 = process_nl_constraints(nlcs, x0)
+    assert all(nlconstr0 == [x0, x0, x0])
+    assert len(constr_func(0)) == 3
     assert all(constr_func(0) == [0, 0, 0])
 
 
-@pytest.mark.parametrize("lb", (-np.inf, [-np.inf], np.array([-np.inf])))
-@pytest.mark.parametrize("ub", (0, [0], np.array([0])))
+@pytest.mark.parametrize("lb", (-np.inf, [-np.inf]*2, np.array([-np.inf]*2)))
+@pytest.mark.parametrize("ub", (0, [0]*2, np.array([0]*2)))
 def test_single_nl_constraint(lb, ub):
     nlc = NonlinearConstraint(lambda x: [x, x], lb=lb, ub=ub)
     x0 = 2.1
-    constr_func = process_nl_constraints([nlc])
+    constr_func, nlconstr0 = process_nl_constraints([nlc], x0)
+    assert all(nlconstr0 == [x0, x0])
+    assert len(constr_func(0)) == 2
     assert all(constr_func(x0) == [x0, x0])
+
+@pytest.mark.parametrize("lb", (-np.inf, [-np.inf]))
+@pytest.mark.parametrize("ub", (0.0, [0.0]))
+def test_length_nlc_fun_not_equal_to_length_lb_ub(lb, ub):
+    if lb == -np.inf and ub == 0.0:
+        return  # No exception here
+    nlc = NonlinearConstraint(lambda x: [x, x], lb=lb, ub=ub)
+    x0 = 2.1
+    with pytest.raises(ValueError) as e_info:
+        process_nl_constraints([nlc], x0)
+    if lb == -np.inf:
+        assert e_info.match("The number of elements in the constraint function's output does not match the number of elements in the upper bound.")
+    else:
+        assert e_info.match("The number of elements in the constraint function's output does not match the number of elements in the lower bound.")
+
+
+def test_length_nlc_fun_ne_to_length_ub():
+    # This specific test is for the case that both lb and ub are vectors, and lb has the correct
+    # length but ub does not.
+    nlc = NonlinearConstraint(lambda x: [x, x], lb=[-np.inf, 0], ub=[0])
+    x0 = 2.1
+    with pytest.raises(ValueError) as e_info:
+        process_nl_constraints([nlc], x0)
+    assert e_info.match("The number of elements in the constraint function's output does not match the number of elements in the upper bound.")
+
+
+def test_lb_neg_inf_ub_inf_raises():
+    nlc = NonlinearConstraint(lambda x: [x, x], lb=-np.inf, ub=np.inf)
+    with pytest.raises(ValueError) as e_info:
+        process_nl_constraints([nlc], 0)
+    assert e_info.match("A NonlinearConstraint was provided without specifying lower or upper bounds")
+
+def test_lb_neg_inf_ub_vector_w_inf_raises():
+    nlc = NonlinearConstraint(lambda x: [x, x], lb=-np.inf, ub=[0, np.inf])
+    with pytest.raises(ValueError) as e_info:
+        process_nl_constraints([nlc], 0)
+    assert e_info.match("A NonlinearConstraint was provided without specifying lower or upper bounds")
+
+def test_lb_vector_with_neg_inf_ub_inf_raises():
+    nlc = NonlinearConstraint(lambda x: [x, x], lb=[-np.inf, 0], ub=np.inf)
+    with pytest.raises(ValueError) as e_info:
+        process_nl_constraints([nlc], 0)
+    assert e_info.match("A NonlinearConstraint was provided without specifying lower or upper bounds")
+
+def test_lb_vector_with_neg_inf_ub_vector_w_inf_at_same_index_raises():
+    nlc = NonlinearConstraint(lambda x: [x, x], lb=[-np.inf, 0], ub=[np.inf, 1])
+    with pytest.raises(ValueError) as e_info:
+        process_nl_constraints([nlc], 0)
+    assert e_info.match("A NonlinearConstraint was provided without specifying lower or upper bounds")