r m x p   toggle line displays

j k   next/prev highlighted chunk

0   (zero) top of page

1   (one) first highlighted chunk

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

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.")

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)))

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(

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()

201 if not isinstance(other, AffineExpression):

202 return NotImplemented

204 return Ellipsoid(self._n, self._A, self._c + other)

206 @convert_operands()

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())