Coverage for picos/expressions/set_ellipsoid.py: 69.91%
113 statements
« prev ^ index » next coverage.py v6.5.0, created at 2023-03-26 07:46 +0000
« prev ^ index » next coverage.py v6.5.0, created at 2023-03-26 07:46 +0000
1# ------------------------------------------------------------------------------
2# Copyright (C) 2020 Maximilian Stahlberg
3#
4# This file is part of PICOS.
5#
6# PICOS is free software: you can redistribute it and/or modify it under the
7# terms of the GNU General Public License as published by the Free Software
8# Foundation, either version 3 of the License, or (at your option) any later
9# version.
10#
11# PICOS is distributed in the hope that it will be useful, but WITHOUT ANY
12# WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR
13# A PARTICULAR PURPOSE. See the GNU General Public License for more details.
14#
15# You should have received a copy of the GNU General Public License along with
16# this program. If not, see <http://www.gnu.org/licenses/>.
17# ------------------------------------------------------------------------------
19"""Implements :class:`Ellipsoid`."""
21import operator
22from collections import namedtuple
24from .. import glyphs
25from ..apidoc import api_end, api_start
26from ..caching import cached_property
27from ..constraints import SOCConstraint
28from .data import convert_operands, cvxopt_inverse
29from .exp_affine import AffineExpression, Constant
30from .expression import Expression
31from .set import Set
33_API_START = api_start(globals())
34# -------------------------------
37class Ellipsoid(Set):
38 r"""An affine transformation of the Euclidean unit ball.
40 :Definition:
42 For :math:`n \in \mathbb{Z}_{\geq 1}`, :math:`A \in \mathbb{R}^{n \times n}`
43 invertible and :math:`c \in \mathbb{R}^n`, an instance of this class
44 represents the set
46 .. math::
48 & \{Ax + c \mid \lVert x \rVert_2 \leq 1\} \\
49 =~& \{x \mid \lVert A^{-1}(x - c) \rVert_2 \leq 1\} \\
50 =~& \{x \mid \lVert x - c \rVert_{(A^{-1})^T A^{-1}} \leq 1\}.
52 Unlike most sets, instances of this class offer a limited set of algebraic
53 operations that are generalized from expressions to sets in the natural way.
54 In particular, you can add or substract constant vectors of matching
55 dimension and apply matrix multiplication from the left hand side, both of
56 which will act on the term :math:`Ax + c` in the definition above.
58 :Example:
60 >>> from picos import Ellipsoid, RealVariable
61 >>> Ellipsoid(3) # Three-dimensional Euclidean unit ball.
62 <Centered Unit Ball: {I·x : ‖x‖ ≤ 1}>
63 >>> Ellipsoid(3, range(9)) # Linear transformation of the unit ball.
64 <Centered Ellipsoid: {[3×3]·x : ‖x‖ ≤ 1}>
65 >>> Ellipsoid(3, "2I", 1) # Offset ball of radius two.
66 <Offset Ellipsoid: {2·I·x + [1] : ‖x‖ ≤ 1}>
67 >>> 2*Ellipsoid(3) + 1 # The same using algebraic operations.
68 <Offset Ellipsoid: {2·I·x + [1] : ‖x‖ ≤ 1}>
69 >>> x = RealVariable("x", 3)
70 >>> (2*x + range(3)) << (4*Ellipsoid(3) + 5) # Constraint creation.
71 <4×1 SOC Constraint: ‖(4·I)^(-1)·(2·x + [3×1] - [5])‖ ≤ 1>
73 .. note::
75 Due to significant differences in scope, :class:`Ellipsoid` is not a
76 superclass of :class:`~.set_ball.Ball` even though both classes can
77 represent Euclidean balls around the origin.
78 """
80 def __init__(self, n, A="I", c=0):
81 """Construct an ellipsoid.
83 :param int n:
84 Dimensionality :math:`n` of the ellipsoid.
86 :param A:
87 Invertible linear transformation matrix :math:`A`.
88 :type A:
89 :class:`~.exp_affine.AffineExpression` or recognized by
90 :func:`~.data.load_data`
92 :param c:
93 The center :math:`c` of the ellispoid.
94 :type c:
95 :class:`~.exp_affine.AffineExpression` or recognized by
96 :func:`~.data.load_data`
98 .. warning::
100 Invertibility of :math:`A` is not checked on instanciation.
101 If :math:`A` is singular, a :exc:`RuntimeError` is raised once the
102 inverse is needed.
103 """
104 if not isinstance(n, int):
105 raise TypeError("Dimensionality must be an integer.")
107 if n < 1:
108 raise ValueError("Dimensionality must be positive.")
110 # Load A.
111 if not isinstance(A, Expression):
112 A = AffineExpression.from_constant(A, shape=(n, n))
113 else:
114 A = A.refined
116 if not isinstance(A, AffineExpression) or not A.constant:
117 raise TypeError("A must be a constant real matrix.")
119 if A.shape != (n, n):
120 raise TypeError("A must be a {} matrix, got {} instead."
121 .format(glyphs.shape((n, n)), glyphs.shape(A.shape)))
123 # Load c.
124 if not isinstance(c, Expression):
125 c = AffineExpression.from_constant(c, shape=n)
126 else:
127 c = c.refined
129 if not isinstance(c, AffineExpression) or not c.constant:
130 raise TypeError("c must be a constant real vector.")
132 if c.shape != (n, 1):
133 raise TypeError(
134 "c must be a {}-dimensional column vector.".format(n))
136 self._n = n
137 self._A = A
138 self._c = c
140 typeStr = (
141 ("Centered " if c.is0 else "Offset ")
142 + ("Unit Ball" if A.isI else "Ellipsoid"))
144 varName = glyphs.free_var_name(A.string + c.string)
145 symbStr = glyphs.set(glyphs.sep(
146 glyphs.clever_add(glyphs.clever_mul(A.string, varName), c.string),
147 glyphs.le(glyphs.norm(varName), 1)))
149 Set.__init__(self, typeStr, symbStr)
151 # --------------------------------------------------------------------------
152 # Properties.
153 # --------------------------------------------------------------------------
155 @property
156 def dim(self):
157 """The dimensionality :math:`n`."""
158 return self._n
160 @property
161 def A(self):
162 """The linear operator matrix :math:`A`."""
163 return self._A
165 @cached_property
166 def Ainv(self):
167 """The inverse linear operator matrix :math:`A^{-1}`."""
168 try:
169 inverse = cvxopt_inverse(self._A.value_as_matrix)
170 except ValueError as error:
171 raise RuntimeError("The matrix A is not invertible.") from error
173 return Constant(glyphs.inverse(self._A.string), inverse, self._A.shape)
175 @property
176 def c(self):
177 """The center point :math:`c`."""
178 return self._c
180 # --------------------------------------------------------------------------
181 # Abstract method implementations.
182 # --------------------------------------------------------------------------
184 def _get_mutables(self):
185 return frozenset()
187 def _replace_mutables(self, mapping):
188 return self
190 Subtype = namedtuple("Subtype", ("dim",))
192 def _get_subtype(self):
193 return self.Subtype(self._n)
195 # --------------------------------------------------------------------------
196 # Algebraic operations.
197 # --------------------------------------------------------------------------
199 @convert_operands()
200 def __add__(self, other):
201 if not isinstance(other, AffineExpression):
202 return NotImplemented
204 return Ellipsoid(self._n, self._A, self._c + other)
206 @convert_operands()
207 def __radd__(self, other):
208 if not isinstance(other, AffineExpression):
209 return NotImplemented
211 return Ellipsoid(self._n, self._A, other + self._c)
213 @convert_operands()
214 def __mul__(self, other):
215 if not isinstance(other, AffineExpression):
216 return NotImplemented
218 if not other.scalar:
219 raise TypeError(
220 "Can only multiply an Ellipsoid from the right with a scalar.")
222 return Ellipsoid(self._n, self._A*other, self._c*other)
224 @convert_operands()
225 def __rmul__(self, other):
226 if not isinstance(other, AffineExpression):
227 return NotImplemented
229 if other.shape not in (self._A.shape, (1, 1)):
230 raise TypeError("Can only multiply a {}-dimensional Ellipsoid from "
231 "the left with a scalar or a {} matrix.".format(
232 self._n, glyphs.shape(self._A.shape)))
234 return Ellipsoid(self._n, other*self._A, other*self._c)
236 @convert_operands()
237 def __truediv__(self, other):
238 if not isinstance(other, AffineExpression):
239 return NotImplemented
241 if not other.scalar:
242 raise TypeError("You may only divide an Ellipsoid by a scalar.")
244 if not other.constant:
245 raise TypeError("You may only divide an Ellipsoid by a constant.")
247 if other.is0:
248 raise ZeroDivisionError("Tried to divide an Ellipsoid by zero.")
250 return Ellipsoid(self._n, self._A/other, self._c/other)
252 # --------------------------------------------------------------------------
253 # Constraint-creating operations.
254 # --------------------------------------------------------------------------
256 @classmethod
257 def _predict(cls, subtype, relation, other):
258 assert isinstance(subtype, cls.Subtype)
260 if relation == operator.__rshift__:
261 if issubclass(other.clstype, AffineExpression):
262 if subtype.dim == other.subtype.dim:
263 return SOCConstraint.make_type(subtype.dim)
265 return NotImplemented
267 def _rshift_implementation(self, element):
268 if isinstance(element, AffineExpression):
269 if len(element) != self._n:
270 raise TypeError("Cannot constrain the {}-dimensional "
271 "expression {} into a {}-dimensional ellipsoid."
272 .format(len(element), element.string, self._n))
274 one = AffineExpression.from_constant(1)
275 return SOCConstraint(self.Ainv*(element.vec - self.c), one)
276 else:
277 return NotImplemented
280# --------------------------------------
281__all__ = api_end(_API_START, globals())