r""" LinQuadDomain definition

Philosophically, LinQuad domains consist of those convex domains specified by a combination
of linear inequality, linear equality, and quadratic inequality constraints.

From a code perspective, this is where a dependency on CVXPY is introduced to
handle domain properities from within this domain

"""

from __future__ import division

import numpy as np
import cvxpy as cp
from .domain import TOL, DEFAULT_CVXPY_KWARGS
from .euclidean import EuclideanDomain
from .tensor import TensorProductDomain

from ..misc import merge

r"""A domain specified by a combination of linear (in)equality constraints and convex quadratic constraints

Here we define a domain that is specified in terms of bound constraints,
linear inequality constraints, linear equality constraints, and quadratic constraints.

.. math::

\mathcal{D} := \left \lbrace
\mathbf{x} : \text{lb} \le \mathbf{x} \le \text{ub}, \
\mathbf{A} \mathbf{x} \le \mathbf{b}, \
\mathbf{A}_{\text{eq}} \mathbf{x} = \mathbf{b}_{\text{eq}}, \
\| \mathbf{L}_i (\mathbf{x} - \mathbf{y}_i)\|_2 \le \rho_i
\right\rbrace \subset \mathbb{R}^m

Parameters
----------
A: array-like (m,n)
Matrix in left-hand side of inequality constraint
b: array-like (m,)
Vector in right-hand side of the ineqaluty constraint
A_eq: array-like (p,n)
Matrix in left-hand side of equality constraint
b_eq: array-like (p,)
Vector in right-hand side of equality constraint
lb: array-like (n,)
Vector of lower bounds
ub: array-like (n,)
Vector of upper bounds
Ls: list of array-likes (p,m)
List of matrices with m columns defining the quadratic constraints
ys: list of array-likes (m,)
rhos: list of positive floats
names: list of strings, optional
Names for each of the parameters in the space
kwargs: dict, optional
Additional parameters to be passed to cvxpy Problem.solve()
"""
def __init__(self, A = None, b = None,
lb = None, ub = None,
A_eq = None, b_eq = None,
Ls = None, ys = None, rhos = None,
names = None, **kwargs):

self.tol = 1e-6
# Determine dimension of space
self._init_dim(lb = lb, ub = ub, A = A, A_eq = A_eq, Ls = Ls)

# Start setting default values
self._lb = self._init_lb(lb)
self._ub = self._init_ub(ub)
self._A, self._b = self._init_ineq(A, b)
self._A_eq, self._b_eq = self._init_eq(A_eq, b_eq)
self._Ls, self._ys, self._rhos = self._init_quad(Ls, ys, rhos)

self._init_names(names)

self.kwargs = merge(DEFAULT_CVXPY_KWARGS, kwargs)

def __str__(self):
ret = "<%s on R^%d" % (self.__class__.__name__, len(self))
if len(self._Ls) > 0:
ret += "; %d quadratic constraints" % (len(self._Ls),)
if self._A.shape > 0:
ret += "; %d linear inequality constraints" % (self._A.shape, )
if self._A_eq.shape > 0:
ret += "; %d linear equality constraints" % (self._A_eq.shape, )
ret +=">"

return ret

################################################################################
# Initialization helpers
################################################################################
def _init_dim(self, lb = None, ub = None, A = None, A_eq = None, Ls = None):
"""determine the dimension of the space we are working on"""
if lb is not None:
if isinstance(lb, (int, float)):
m = 1
else:
m = len(lb)
elif ub is not None:
if isinstance(ub, (int, float)):
m = 1
else:
m = len(ub)
elif A is not None:
m = len(A)
elif A_eq is not None:
m = len(A_eq)
elif Ls is not None:
m = len(Ls)
else:
raise Exception("Could not determine dimension of space")

self._dimension = m

def _init_lb(self, lb):
if lb is None:
return -np.inf*np.ones(len(self))
else:
if isinstance(lb, (int, float)):
lb = [lb]
assert len(lb) == len(self), "Lower bound has wrong dimensions"
return np.array(lb)

def _init_ub(self, ub):
if ub is None:
return np.inf*np.ones(len(self))
else:
if isinstance(ub, (int, float)):
ub = [ub]
assert len(ub) == len(self), "Upper bound has wrong dimensions"
return np.array(ub)

def _init_ineq(self, A, b):
if A is None and b is None:
A = np.zeros((0,len(self)))
b = np.zeros((0,))
elif A is not None and b is not None:
A = np.array(A)
b = np.array(b)
if len(b.shape) == 0:
b = b.reshape(1)

assert len(b.shape) == 1, "b must have only one dimension"

if len(A.shape) == 1 and len(b) == 1:
A = A.reshape(1,-1)

assert A.shape == len(self), "A has wrong number of columns"
assert A.shape == b.shape, "The number of rows of A and b do not match"
else:
raise AssertionError("If using inequality constraints, both A and b must be specified")
return A, b

def _init_eq(self, A_eq, b_eq):
if A_eq is None and b_eq is None:
A_eq = np.zeros((0,len(self)))
b_eq = np.zeros((0,))
elif A_eq is not None and b_eq is not None:
A_eq = np.array(A_eq)
b_eq = np.array(b_eq)
if len(b_eq.shape) == 0:
b_eq = b_eq.reshape(1)

assert len(b_eq.shape) == 1, "b_eq must have only one dimension"

if len(A_eq.shape) == 1 and len(b_eq) == 1:
A_eq = A_eq.reshape(1,-1)

assert A_eq.shape == len(self), "A_eq has wrong number of columns"
assert A_eq.shape == b_eq.shape, "The number of rows of A_eq and b_eq do not match"
else:
raise AssertionError("If using equality constraints, both A_eq and b_eq must be specified")

return A_eq, b_eq

if Ls is None and ys is None and rhos is None:
_Ls = []
_ys = []
_rhos = []
elif Ls is not None and ys is not None and rhos is not None:
assert len(Ls) == len(ys) == len(rhos), "Length of all quadratic constraints must be the same"

_Ls = []
_ys = []
_rhos = []
for L, y, rho in zip(Ls, ys, rhos):
assert len(L) == len(self), "dimension of L doesn't match the domain"
assert len(y) == len(self), "Dimension of center doesn't match the domain"
assert rho > 0, "Radius must be positive"
_Ls.append(np.array(L))
_ys.append(np.array(y))
_rhos.append(rho)
# TODO: If constraint is rank-1, should we implicitly convert to a linear inequality constriant
else:
raise AssertionError("If providing quadratic constraint, each of Ls, ys, and rhos must be defined")
return _Ls, _ys, _rhos

################################################################################
# Simple properties
################################################################################
def __len__(self): return self._dimension

@property
def lb(self): return self._lb

@property
def ub(self): return self._ub

@property
def A(self): return self._A

@property
def b(self): return self._b

@property
def A_eq(self): return self._A_eq

@property
def b_eq(self): return self._b_eq

@property
def Ls(self): return self._Ls

@property
def ys(self): return self._ys

@property
def rhos(self): return self._rhos

################################################################################
# Normalization
################################################################################
def _normalized_domain(self, **kwargs):
names_norm = [name + ' (normalized)' for name in self.names]

return LinQuadDomain(lb = self.lb_norm, ub = self.ub_norm, A = self.A_norm, b = self.b_norm,
A_eq = self.A_eq_norm, b_eq = self.b_eq_norm, Ls = self.Ls_norm, ys = self.ys_norm, rhos = self.rhos_norm,
names = names_norm, **merge(self.kwargs, kwargs))

def _isinside(self, X, tol = TOL):
return self._isinside_bounds(X, tol = tol) & self._isinside_ineq(X, tol = tol) & self._isinside_eq(X, tol = tol) & self._isinside_quad(X, tol = tol)

def _extent(self, x, p):
# Check that direction satisfies equality constraints to a tolerance
if self.A_eq.shape == 0 or np.all(np.abs(self.A_eq.dot(p) ) < self.tol):
return min(self._extent_bounds(x, p), self._extent_ineq(x, p), self._extent_quad(x, p))
else:
return 0.

################################################################################
# Convex Solver Functions
################################################################################

def _build_constraints_norm(self, x_norm):
r""" Build the constraints corresponding to the domain given a vector x
"""
constraints = []

# Numerical issues emerge with unbounded constraints
I = np.isfinite(self.lb_norm)
if np.sum(I) > 0:
constraints.append( self.lb_norm[I] <= x_norm[I])

I = np.isfinite(self.ub_norm)
if np.sum(I) > 0:
constraints.append( x_norm[I] <= self.ub_norm[I])

if self.A.shape > 0:
constraints.append( x_norm.__rmatmul__(self.A_norm) <= self.b_norm)
if self.A_eq.shape > 0:
constraints.append( x_norm.__rmatmul__(self.A_eq_norm) == self.b_eq_norm)

for L, y, rho in zip(self.Ls_norm, self.ys_norm, self.rhos_norm):
if len(L) > 1:
constraints.append( cp.norm(x_norm.__rmatmul__(L) - L.dot(y)) <= rho )
elif len(L) == 1:
constraints.append( cp.norm(L*x_norm - L.dot(y)) <= rho)

return constraints

def _build_constraints(self, x):
r""" Build the constraints corresponding to the domain given a vector x
"""
constraints = []

# Numerical issues emerge with unbounded constraints
I = np.isfinite(self.lb)
if np.sum(I) > 0:
constraints.append( self.lb[I] <= x[I])

I = np.isfinite(self.ub)
if np.sum(I) > 0:
constraints.append( x[I] <= self.ub[I])

if self.A.shape > 0:
constraints.append( x.__rmatmul__(self.A) <= self.b)
if self.A_eq.shape > 0:
constraints.append( x.__rmatmul__(self.A_eq) == self.b_eq)

for L, y, rho in zip(self.Ls, self.ys, self.rhos):
if len(L) > 1:
constraints.append( cp.norm(x.__rmatmul__(L) - L.dot(y)) <= rho )
elif len(L) == 1:
constraints.append( cp.norm(L*x - L.dot(y)) <= rho)

return constraints

################################################################################
#
################################################################################

def add_constraints(self, A = None, b = None, lb = None, ub = None, A_eq = None, b_eq = None,
Ls = None, ys = None, rhos = None):
r"""Add new constraints to the domain
"""
lb = self._init_lb(lb)
ub = self._init_ub(ub)
A, b = self._init_ineq(A, b)
A_eq, b_eq = self._init_eq(A_eq, b_eq)
Ls, ys, rhos = self._init_quad(Ls, ys, rhos)

# Update constraints
lb = np.maximum(lb, self.lb)
ub = np.minimum(ub, self.ub)

A = np.vstack([self.A, A])
b = np.hstack([self.b, b])

A_eq = np.vstack([self.A_eq, A_eq])
b_eq = np.hstack([self.b_eq, b_eq])

Ls = self.Ls + Ls
ys = self.ys + ys
rhos = self.rhos + rhos

if len(Ls) > 0:
return LinQuadDomain(lb = lb, ub = ub, A = A, b = b, A_eq = A_eq, b_eq = b_eq,
Ls = Ls, ys = ys, rhos = rhos, names = self.names)
elif len(b) > 0 or len(b_eq) > 0:
from .linineq import LinIneqDomain
return LinIneqDomain(lb = lb, ub = ub, A = A, b = b, A_eq = A_eq, b_eq = b_eq, names = self.names)
else:
from .box import BoxDomain
return BoxDomain(lb = lb, ub = ub, names = self.names)

def __and__(self, other):