from __future__ import division, print_function, absolute_import
import numpy as np
import scipy.sparse as spc
from scipy.sparse.linalg import LinearOperator
from ._numdiff import approx_derivative
from warnings import warn
__all__ = ['NonlinearConstraint',
'LinearConstraint',
'BoxConstraint']
class NonlinearConstraint:
"""Nonlinear constraint
Parameters
----------
fun : callable
The function defining the constraint.
fun(x) -> array_like, shape (m,)
where x is a (n,) ndarray and ``m``
is the number of constraints.
kind : {str, tuple}
Specifies the type of contraint. Options for this
parameter are:
- ('interval', lb, ub) for a constraint of the type:
lb <= fun(x) <= ub
- ('greater', lb) for a constraint of the type:
fun(x) >= lb
- ('less', ub) for a constraint of the type:
fun(x) <= ub
- ('equals', c) for a constraint of the type:
fun(x) == c
- ('greater',) for a constraint of the type:
fun(x) >= 0
- ('less',) for a constraint of the type:
fun(x) <= 0
- ('equals',) for a constraint of the type:
fun(x) == 0
where ``lb``, ``ub`` and ``c`` are (m,) ndarrays or
scalar values. In the latter case, the same value
will be repeated for all the constraints.
jac : callable
Jacobian Matrix:
jac(x) -> {ndarray, sparse matrix}, shape (m, n)
where x is a (n,) ndarray.
hess : {callable, '2-point', '3-point', 'cs', None}
Method for computing the Hessian matrix. The keywords
select a finite difference scheme for numerical
estimation. The scheme '3-point' is more accurate, but requires
twice as much operations compared to '2-point' (default). The
scheme 'cs' uses complex steps, and while potentially the most
accurate, it is applicable only when `fun` correctly handles
complex inputs and can be analytically continued to the complex
plane. If it is a callable, it should return the
Hessian matrix of `dot(fun, v)`:
hess(x, v) -> {LinearOperator, sparse matrix, ndarray}, shape (n, n)
where x is a (n,) ndarray and v is a (m,) ndarray. When ``hess``
is None it considers the hessian is an matrix filled with zeros.
enforce_feasibility : {list of boolean, boolean}, optional
Specify if the constraint must be feasible along the iterations.
If ``True`` all the iterates generated by the optimization
algorithm need to be feasible in respect to a constraint. If ``False``
this is not needed. A list can be passed to specify element-wise
each constraints needs to stay feasible along the iterations and
each does not. Alternatively, a single boolean can be used to
specify the feasibility required of all constraints. By default it
is False.
"""
def __init__(self, fun, kind, jac, hess='2-point', enforce_feasibility=False):
self._fun = fun
self.kind = kind
self._jac = jac
self._hess = hess
self.enforce_feasibility = enforce_feasibility
self.isinitialized = False
def evaluate_and_initialize(self, x0, sparse_jacobian=None):
x0 = np.atleast_1d(x0).astype(float)
f0 = np.atleast_1d(self._fun(x0))
v0 = np.zeros_like(f0)
J0 = self._jac(x0)
def fun_wrapped(x):
return np.atleast_1d(self._fun(x))
if sparse_jacobian or (sparse_jacobian is None and spc.issparse(J0)):
def jac_wrapped(x):
return spc.csr_matrix(self._jac(x))
self.sparse_jacobian = True
self.J0 = spc.csr_matrix(J0)
else:
def jac_wrapped(x):
J = self._jac(x)
if spc.issparse(J):
return J.toarray()
else:
return np.atleast_2d(J)
self.sparse_jacobian = False
if spc.issparse(J0):
self.J0 = J0.toarray()
else:
self.J0 = np.atleast_2d(J0)
if callable(self._hess):
H0 = self._hess(x0, v0)
if spc.issparse(H0):
H0 = spc.csr_matrix(H0)
def hess_wrapped(x, v):
return spc.csr_matrix(self._hess(x, v))
elif isinstance(H0, LinearOperator):
def hess_wrapped(x, v):
return self._hess(x, v)
else:
H0 = np.atleast_2d(np.asarray(H0))
def hess_wrapped(x, v):
return np.atleast_2d(np.asarray(self._hess(x, v)))
elif self._hess in ('2-point', '3-point', 'cs'):
approx_method = self._hess
def jac_dot_v(x, v):
J = jac_wrapped(x)
return J.T.dot(v)
def hess_wrapped(x, v):
return approx_derivative(jac_dot_v, x, approx_method,
as_linear_operator=True,
args=(v,))
else:
hess_wrapped = self._hess
self.fun = fun_wrapped
self.jac = jac_wrapped
self.hess = hess_wrapped
self.x0 = x0
self.f0 = f0
self.n = x0.size
self.m = f0.size
self.kind = _check_kind(self.kind, self.m)
self.enforce_feasibility \
= _check_enforce_feasibility(self.enforce_feasibility, self.m)
if not _is_feasible(self.kind, self.enforce_feasibility, f0):
raise ValueError("Unfeasible initial point. "
"Either set ``enforce_feasibility=False`` or "
"choose a new feasible initial point ``x0``.")
self.isinitialized = True
return x0
class LinearConstraint:
"""Linear constraint.
Parameters
----------
A : {ndarray, sparse matrix}, shape (m, n)
Matrix for the linear constraint.
kind : {str, tuple}
Specifies the type of contraint. Options for this
parameter are:
- ('interval', lb, ub) for a constraint of the type:
lb <= A x <= ub
- ('greater', lb) for a constraint of the type:
A x >= lb
- ('less', ub) for a constraint of the type:
A x <= ub
- ('equals', c) for a constraint of the type:
A x == c
- ('greater',) for a constraint of the type:
A x >= 0
- ('less',) for a constraint of the type:
A x <= 0
- ('equals',) for a constraint of the type:
A x == 0
where ``lb``, ``ub`` and ``c`` are (m,) ndarrays or
scalar values. In the latter case, the same value
will be repeated for all the constraints.
enforce_feasibility : {list of boolean, boolean}, optional
Specify if the constraint must be feasible along the iterations.
If ``True`` all the iterates generated by the optimization
algorithm need to be feasible in respect to a constraint. If ``False``
this is not needed. A list can be passed to specify element-wise
each constraints needs to stay feasible along the iterations and
each does not. Alternatively, a single boolean can be used to
specify the feasibility required of all constraints. By default it
is False.
"""
def __init__(self, A, kind, enforce_feasibility=False):
self.A = A
self.kind = kind
self.enforce_feasibility = enforce_feasibility
self.isinitialized = False
def evaluate_and_initialize(self, x0, sparse_jacobian=None):
if sparse_jacobian or (sparse_jacobian is None
and spc.issparse(self.A)):
self.A = spc.csr_matrix(self.A)
self.sparse_jacobian = True
else:
if spc.issparse(self.A):
self.A = self.A.toarray()
else:
self.A = np.atleast_2d(self.A)
self.sparse_jacobian = False
x0 = np.atleast_1d(x0).astype(float)
f0 = self.A.dot(x0)
J0 = self.A
self.x0 = x0
self.f0 = f0
self.J0 = J0
self.n = x0.size
self.m = f0.size
self.kind = _check_kind(self.kind, self.m)
self.enforce_feasibility \
= _check_enforce_feasibility(self.enforce_feasibility, self.m)
if not _is_feasible(self.kind, self.enforce_feasibility, f0):
raise ValueError("Unfeasible initial point. "
"Either set ``enforce_feasibility=False`` or "
"choose a new feasible initial point ``x0``.")
self.isinitialized = True
return x0
def to_nonlinear(self):
if not self.isinitialized:
raise RuntimeError("Trying to convert uninitialized constraint.")
def fun(x):
return self.A.dot(x)
def jac(x):
return self.A
# Build Constraints
nonlinear = NonlinearConstraint(fun, self.kind, jac, None,
self.enforce_feasibility)
nonlinear.isinitialized = True
nonlinear.m = self.m
nonlinear.n = self.n
nonlinear.sparse_jacobian = self.sparse_jacobian
nonlinear.fun = fun
nonlinear.jac = jac
nonlinear.hess = None
nonlinear.x0 = self.x0
nonlinear.f0 = self.f0
nonlinear.J0 = self.J0
return nonlinear
class BoxConstraint:
"""Box constraint.
Parameters
----------
kind : tuple
Specifies the type of contraint. Options for this
parameter are:
- ('interval', lb, ub) for a constraint of the type:
lb <= A x <= ub
- ('greater', lb) for a constraint of the type:
A x >= lb
- ('less', ub) for a constraint of the type:
A x <= ub
- ('equals', c) for a constraint of the type:
A x == c
- ('greater',) for a constraint of the type:
A x >= 0
- ('less',) for a constraint of the type:
A x <= 0
- ('equals',) for a constraint of the type:
A x == 0
where ``lb``, ``ub`` and ``c`` are (m,) ndarrays or
scalar values. In the latter case, the same value
will be repeated for all the constraints.
enforce_feasibility : {list of boolean, boolean}, optional
Specify if the constraint must be feasible along the iterations.
If ``True`` all the iterates generated by the optimization
algorithm need to be feasible in respect to a constraint. If ``False``
this is not needed. A list can be passed to specify element-wise
each constraints needs to stay feasible along the iterations and
each does not. Alternatively, a single boolean can be used to
specify the feasibility required of all constraints. By default it
is False.
"""
def __init__(self, kind, enforce_feasibility=False):
self.kind = kind
self.enforce_feasibility = enforce_feasibility
self.isinitialized = False
def evaluate_and_initialize(self, x0, sparse_jacobian=None):
x0 = np.atleast_1d(x0).astype(float)
f0 = x0
self.n = x0.size
self.m = f0.size
if sparse_jacobian or sparse_jacobian is None:
J0 = spc.eye(self.n).tocsr()
self.sparse_jacobian = True
else:
J0 = np.eye(self.n)
self.sparse_jacobian = False
self.J0 = J0
self.kind = _check_kind(self.kind, self.m)
self.enforce_feasibility \
= _check_enforce_feasibility(self.enforce_feasibility, self.m)
self.isinitialized = True
if not _is_feasible(self.kind, self.enforce_feasibility, f0):
warn("The initial point was changed in order "
"to stay inside box constraints.")
x0_new = _reinforce_box_constraint(self.kind,
self.enforce_feasibility,
x0)
self.x0 = x0_new
self.f0 = x0_new
return x0_new
else:
self.x0 = x0
self.f0 = f0
return x0
def to_linear(self):
if not self.isinitialized:
raise RuntimeError("Trying to convert uninitialized constraint.")
# Build Constraints
linear = LinearConstraint(self.J0, self.kind,
self.enforce_feasibility)
linear.isinitialized = True
linear.m = self.m
linear.n = self.n
linear.sparse_jacobian = self.sparse_jacobian
linear.x0 = self.x0
linear.f0 = self.f0
linear.J0 = self.J0
return linear
def to_nonlinear(self):
if not self.isinitialized:
raise RuntimeError("Trying to convert uninitialized constraint.")
return self.to_linear().to_nonlinear()
# ************************************************************ #
# ********** Auxiliar Functions ********** #
# ************************************************************ #
def _check_kind(kind, m):
if not isinstance(kind, (tuple, list, str)):
raise ValueError("The parameter `kind` should be a tuple, "
" a list, or a string.")
if isinstance(kind, str):
kind = (kind,)
if len(kind) == 0:
raise ValueError("The parameter `kind` should not be empty.")
n_args = len(kind)
keyword = kind[0]
if keyword not in ("greater", "less", "equals", "interval"):
raise ValueError("Keyword `%s` not available." % keyword)
if n_args in (1, 2) and keyword not in ("greater", "less", "equals") \
or n_args == 3 and keyword not in ("interval"):
raise ValueError("Invalid `kind` format.")
if n_args == 1:
kind = (keyword, 0)
if keyword in ("greater", "less", "equals"):
c = np.asarray(kind[1], dtype=float)
if np.size(c) not in (1, m):
if keyword == "greater":
raise ValueError("`lb` has the wrong dimension.")
if keyword == "less":
raise ValueError("`ub` has the wrong dimension.")
if keyword == "equals":
raise ValueError("`c` has the wrong dimension.")
c = np.resize(c, m)
return (keyword, c)
elif keyword == "interval":
lb = np.asarray(kind[1], dtype=float)
if np.size(lb) not in (1, m):
raise ValueError("`lb` has the wrong dimension.")
lb = np.resize(lb, m)
ub = np.asarray(kind[2], dtype=float)
if np.size(ub) not in (1, m):
raise ValueError("`ub` has the wrong dimension.")
ub = np.resize(ub, m)
if (lb > ub).any():
raise ValueError("lb[i] > ub[i].")
return (keyword, lb, ub)
def _check_enforce_feasibility(enforce_feasibility, m):
if isinstance(enforce_feasibility, bool):
enforce_feasibility = np.full(m,
enforce_feasibility,
dtype=bool)
else:
enforce_feasibility = np.array(enforce_feasibility,
dtype=bool)
if enforce_feasibility.size != m:
raise ValueError("The parameter 'enforce_feasibility' "
"has the wrong number of elements.")
return enforce_feasibility
def _is_feasible(kind, enforce_feasibility, f0):
keyword = kind[0]
if keyword == "equals":
lb = np.asarray(kind[1], dtype=float)
ub = np.asarray(kind[1], dtype=float)
elif keyword == "greater":
lb = np.asarray(kind[1], dtype=float)
ub = np.full_like(lb, np.inf, dtype=float)
elif keyword == "less":
ub = np.asarray(kind[1], dtype=float)
lb = np.full_like(ub, -np.inf, dtype=float)
elif keyword == "interval":
lb = np.asarray(kind[1], dtype=float)
ub = np.asarray(kind[2], dtype=float)
else:
raise RuntimeError("Never be here.")
return ((lb[enforce_feasibility] <= f0[enforce_feasibility]).all()
and (f0[enforce_feasibility] <= ub[enforce_feasibility]).all())
def _reinforce_box_constraint(kind, enforce_feasibility, x0,
relative_tolerance=0.01,
absolute_tolerance=0.01):
"""Reinforce box constraint"""
x0 = np.copy(np.asarray(x0, dtype=float))
keyword = kind[0]
if keyword == "greater":
lb = np.asarray(kind[1], dtype=float)
ub = np.full_like(lb, np.inf, dtype=float)
elif keyword == "less":
ub = np.asarray(kind[1], dtype=float)
lb = np.full_like(ub, -np.inf, dtype=float)
elif keyword == "interval":
lb = np.asarray(kind[1], dtype=float)
ub = np.asarray(kind[2], dtype=float)
x0_new = np.copy(x0)
for i in range(np.size(x0)):
if enforce_feasibility[i]:
if not np.isinf(lb[i]):
lower_bound = min(lb[i]+absolute_tolerance,
lb[i]+relative_tolerance*(ub[i]-lb[i]))
x0_new[i] = max(x0_new[i], lower_bound)
if not np.isinf(ub[i]):
upper_bound = max(ub[i]-absolute_tolerance,
ub[i]-relative_tolerance*(ub[i]-lb[i]))
x0_new[i] = min(x0_new[i], upper_bound)
return x0_new
