Several things:

-Remove custom LC/NLC
-Project x0 for all problems without NLC
-nan_to_num now keeps +/- inf as +/- inf
-cleanup of _common.py

python/prima/__init__.py

@@ -1,40 +1,34 @@
 from ._prima import minimize as _minimize, __version__, PRIMAMessage
-from ._nonlinear_constraints import NonlinearConstraint, process_nl_constraints
+# 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
 from ._linear_constraints import (
-    LinearConstraint,
     process_single_linear_constraint,
     process_multiple_linear_constraints,
     separate_LC_into_eq_and_ineq,
 )
-from ._bounds import process_bounds, Bounds
+from ._bounds import process_bounds
 from enum import Enum
 import numpy as np
 from ._common import _project
 
 
 class ConstraintType(Enum):
-    LINEAR_NATIVE = 5
-    NONLINEAR_NATIVE = 10
-    LINEAR_NONNATIVE = 15
-    NONLINEAR_NONNATIVE = 20
-    LINEAR_DICT = 25
-    NONLINEAR_DICT = 30
+    LINEAR_OBJECT = 5
+    NONLINEAR_OBJECT = 10
+    LINEAR_DICT = 15
+    NONLINEAR_DICT = 20
 
 
 def get_constraint_type(constraint):
-    # Make sure the test for native is first, since the hasattr tests will also pass for native constraints
-    if isinstance(constraint, LinearConstraint):
-        return ConstraintType.LINEAR_NATIVE
-    elif isinstance(constraint, NonlinearConstraint):
-        return ConstraintType.NONLINEAR_NATIVE
-    elif isinstance(constraint, dict) and ("A" in constraint) and ("lb" in constraint) and ("ub" in 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_NONNATIVE
+        return ConstraintType.LINEAR_OBJECT
     elif hasattr(constraint, "fun") and hasattr(constraint, "lb") and hasattr(constraint, "ub"):
-        return ConstraintType.NONLINEAR_NONNATIVE
+        return ConstraintType.NONLINEAR_OBJECT
     else:
         raise ValueError("Constraint type not recognized")
 
@@ -54,15 +48,9 @@ def process_constraints(constraints, x0, options):
     nonlinear_constraints = []
     for constraint in constraints:
         constraint_type = get_constraint_type(constraint)
-        if constraint_type in (
-            ConstraintType.LINEAR_NATIVE,
-            ConstraintType.LINEAR_NONNATIVE,
-        ):
+        if constraint_type is ConstraintType.LINEAR_OBJECT:
             linear_constraints.append(constraint)
-        elif constraint_type in (
-            ConstraintType.NONLINEAR_NATIVE,
-            ConstraintType.NONLINEAR_NONNATIVE,
-        ):
+        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"]))
@@ -128,15 +116,15 @@ def minimize(fun, x0, args=(), method=None, bounds=None, constraints=(), callbac
 
     lb, ub = process_bounds(bounds, lenx0)
 
+    # Project x0 onto the feasible set
+    if nonlinear_constraint is None:
+        result = _project(x0, lb, ub, {"linear": linear_constraint, "nonlinear": None})
+        x0 = result.x
+
     if linear_constraint is not None:
-        # this function doesn't take nonlinear constraints into account at this time
-        x0 = _project(x0, lb, ub, {"linear": linear_constraint, "nonlinear": None})
         A_eq, b_eq, A_ineq, b_ineq = separate_LC_into_eq_and_ineq(linear_constraint)
     else:
-        A_eq = None
-        b_eq = None
-        A_ineq = None
-        b_ineq = None
+        A_eq, b_eq, A_ineq, b_ineq = None, None, None, None
 
     if nonlinear_constraint is not None:
         # PRIMA prefers -inf < f(x) <= 0, so we need to modify the nonlinear constraint accordingly

python/prima/_bounds.py

@@ -32,8 +32,8 @@ def process_bounds(bounds, lenx0):
 
     # This will handle an object of type scipy.optimize.Bounds or similar
     if hasattr(bounds, 'lb') and hasattr(bounds, 'ub'):
-        lb = np.nan_to_num(np.array(bounds.lb, dtype=np.float64), nan=-np.inf)
-        ub = np.nan_to_num(np.array(bounds.ub, dtype=np.float64), nan=np.inf)
+        lb = np.nan_to_num(np.array(bounds.lb, dtype=np.float64), nan=-np.inf, posinf=np.inf, neginf=-np.inf)
+        ub = np.nan_to_num(np.array(bounds.ub, dtype=np.float64), nan=np.inf, posinf=np.inf, neginf=-np.inf)
         lb = np.concatenate((lb, -np.inf*np.ones(lenx0 - len(lb))))
         ub = np.concatenate((ub, np.inf*np.ones(lenx0 - len(ub))))
         return lb, ub

python/prima/_common.py

@@ -1,14 +1,12 @@
 import numpy as np
 from ._linear_constraints import LinearConstraint
+from scipy.optimize import OptimizeResult
 
 # All the accepted scalar types; np.generic correspond to all NumPy types.
 scalar_types = (int, float, np.generic)
 eps = np.finfo(np.float64).eps
-solver_list = ['uobyqa', 'newuoa', 'bobyqa', 'lincoa', 'cobyla']
-invoker_list = solver_list[:]
-# invoker_list.append('pdfo')
 
-def _project(x0, lb, ub, constraints, options=None):
+def _project(x0, lb, ub, constraints):
     """Projection of the initial guess onto the feasible set.
 
     Parameters
@@ -27,11 +25,10 @@ def _project(x0, lb, ub, constraints, options=None):
             nonlinear: dict
                 The nonlinear constraints of the problem. When ``_project`` is called, the nonlinear constraints are
                 None.
-    options: dict, optional
 
     Returns
     -------
-    result: x0
+    result: OptimizeResult
         The result of the projection.
 
     Authors
@@ -43,13 +40,7 @@ def _project(x0, lb, ub, constraints, options=None):
 
     Dedicated to the late Professor M. J. D. Powell FRS (1936--2015).
     """
-    # possible solvers
-    # fun_name = stack()[0][3]  # name of the current function
-    # local_invoker_list = ['prepdfo']
-    # if len(stack()) < 3 or stack()[1][3].lower() not in local_invoker_list:
-    #     raise SystemError('`{}` should only be called by {}'.format(fun_name, ', '.join(invoker_list)))
-    # invoker = stack()[1][3].lower()
-    invoker = ''
+    invoker = 'prima'
 
     # Validate x0.
     if isinstance(x0, scalar_types):
@@ -100,17 +91,13 @@ def _project(x0, lb, ub, constraints, options=None):
         # the nonlinear constraints will not be taken into account in this function and are, therefore, not validated
         raise ValueError('{}: UNEXPECTED ERROR: The constraints are ill-defined.'.format(invoker))
 
-    # Validate options
-    if options is not None and not isinstance(options, dict):
-        raise ValueError('{}: UNEXPECTED ERROR: The options should be a dictionary.'.format(invoker))
-
     max_con = 1e20  # Decide whether an inequality constraint can be ignored
 
     # Project onto the feasible set.
     if constraints['linear'] is None:
         # Direct projection onto the bound constraints
         x_proj = np.nanmin((np.nanmax((x0_c, lb_c), axis=0), ub_c), axis=0)
-        return x_proj
+        return OptimizeResult(x=x_proj)
     elif all(np.less_equal(np.abs(constraints['linear'].ub - constraints['linear'].lb), eps)) and \
             np.max(lb_c) <= -max_con and np.min(ub_c) >= max_con:
         # The linear constraints are all equality constraints. The projection can therefore be done by solving the
@@ -123,12 +110,10 @@ def _project(x0, lb, ub, constraints, options=None):
         # than max_con, they will be reduced to this bound.
         x_proj = np.nanmin((np.nanmax((x0_c + xi, lb_c), axis=0), ub_c), axis=0)
 
-        return x_proj
+        return OptimizeResult(x=x_proj)
 
     if constraints['linear'] is not None:
         try:
-            # TODO: Ideally we'd like to not depend on scipy in the prima package, so in the future we might want to bring in
-            # SLSQP ourselves
             # Project the initial guess onto the linear constraints via SciPy.
             from scipy.optimize import minimize
             from scipy.optimize import Bounds as ScipyBounds
@@ -167,17 +152,16 @@ def _project(x0, lb, ub, constraints, options=None):
             # Perform the actual projection.
             ax_ineq = np.dot(pc_args_ineq['A'], x0_c)
             ax_eq = np.dot(pc_args_eq['A'], x0_c)
-            if np.all(pc_args_ineq['lb'] <= ax_ineq) and np.all(ax_ineq <= pc_args_ineq['ub']) and \
-                    np.all(ax_eq == pc_args_eq['lb']) and \
-                    np.all(lb_c <= x0_c) and np.all(x0_c <= ub_c):
-                # Do not perform any projection if the initial guess is feasible.
-                return x0_c
-            else:
-                res = minimize(lambda x: np.dot(x - x0_c, x - x0_c) / 2, x0_c, jac=lambda x: (x - x0_c),
+            if np.greater(ax_ineq, pc_args_ineq['ub']).any() or np.greater(pc_args_ineq['lb'], ax_ineq).any() or \
+                    np.not_equal(ax_eq, pc_args_eq['lb']).any() or \
+                    np.greater(x0_c, ub_c).any() or np.greater(lb_c, x0_c).any():
+                return minimize(lambda x: np.dot(x - x0_c, x - x0_c) / 2, x0_c, jac=lambda x: (x - x0_c),
                                 bounds=ScipyBounds(lb_c, ub_c), constraints=project_constraints)
-                return res.x
+            else:
+                # Do not perform any projection if the initial guess is feasible.
+                return OptimizeResult(x=x0_c)
 
         except ImportError:
-            return x0_c
+            return OptimizeResult(x=x0_c)
 
-    return x0_c
+    return OptimizeResult(x=x0_c)

python/prima/_linear_constraints.py

@@ -1,38 +1,5 @@
 import numpy as np
-
-
-class LinearConstraint:
-  # Defaults for lb/ub are -inf and inf so that only one needs to be provided and the
-  # other will have no impact
-  def __init__(self, A, lb=-np.inf, ub=np.inf):
-    self.A = A
-    self.lb = lb
-    self.ub = ub
-    # A must be a scalar, a list, or a numpy array
-    # If it's a list, we assume that it's basically a 1-row matrix
-    # If it's a numpy array with only 1 dimension we force it to be a 1-row matrix
-    if isinstance(self.A, int) or isinstance(self.A, float):
-      self.A = np.array([[self.A]])
-    elif isinstance(self.A, list) or isinstance(self.A, np.ndarray):
-      self.A = np.atleast_2d(self.A)
-    else:
-      raise("A must be a scalar, list, or numpy array")
-
-    num_constraints = self.A.shape[0]
-
-    # bounds must be a scalar or a 1D list/array of length equal to number of rows of A
-    def process_bound(bound, name):
-      if isinstance(bound, int) or isinstance(bound, float):
-        bound = np.array([bound]*num_constraints)
-      elif isinstance(bound, list) or isinstance(bound, np.ndarray):
-        bound = np.array(bound)
-        assert bound.ndim == 1 and len(bound) == num_constraints, f"{name} must be a scalar or a 1D list/array of length equal to the number of rows of A"
-      else:
-        raise(f"{name} must be a scalar, list, or numpy array")
-      return bound
-
-    self.lb = process_bound(self.lb, 'lb')
-    self.ub = process_bound(self.ub, 'ub')
+from scipy.optimize import LinearConstraint
 
 def process_single_linear_constraint(constraint):
   # Convert lb and ub to vectors if they are scalars

python/prima/_nonlinear_constraints.py

@@ -1,13 +1,5 @@
 import numpy as np
-
-from warnings import warn  # TODO: Need to determine the text of the warning and actually warn about it
-
-
-class NonlinearConstraint(object):
-    def __init__(self, fun, lb=-np.inf, ub=0):
-        self.fun = fun
-        self.lb = lb
-        self.ub = ub
+from scipy.optimize import NonlinearConstraint
 
 
 def process_nl_constraints(x0, nlcs, options):