Hide keyboard shortcuts

Hot-keys on this page

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 

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

18 

19"""Implements :class:`Ellipsoid`.""" 

20 

21import operator 

22from collections import namedtuple 

23 

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 

32 

33_API_START = api_start(globals()) 

34# ------------------------------- 

35 

36 

37class Ellipsoid(Set): 

38 r"""An affine transformation of the Euclidean unit ball. 

39 

40 :Definition: 

41 

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 

45 

46 .. math:: 

47 

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\}. 

51 

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. 

57 

58 :Example: 

59 

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> 

72 

73 .. note:: 

74 

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

79 

80 def __init__(self, n, A="I", c=0): 

81 """Construct an ellipsoid. 

82 

83 :param int n: 

84 Dimensionality :math:`n` of the ellipsoid. 

85 

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` 

91 

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` 

97 

98 .. warning:: 

99 

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

106 

107 if n < 1: 

108 raise ValueError("Dimensionality must be positive.") 

109 

110 # Load A. 

111 if not isinstance(A, Expression): 

112 A = AffineExpression.from_constant(A, shape=(n, n)) 

113 else: 

114 A = A.refined 

115 

116 if not isinstance(A, AffineExpression) or not A.constant: 

117 raise TypeError("A must be a constant real matrix.") 

118 

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

122 

123 # Load c. 

124 if not isinstance(c, Expression): 

125 c = AffineExpression.from_constant(c, shape=n) 

126 else: 

127 c = c.refined 

128 

129 if not isinstance(c, AffineExpression) or not c.constant: 

130 raise TypeError("c must be a constant real vector.") 

131 

132 if c.shape != (n, 1): 

133 raise TypeError( 

134 "c must be a {}-dimensional column vector.".format(n)) 

135 

136 self._n = n 

137 self._A = A 

138 self._c = c 

139 

140 typeStr = ( 

141 ("Centered " if c.is0 else "Offset ") 

142 + ("Unit Ball" if A.isI else "Ellipsoid")) 

143 

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

148 

149 Set.__init__(self, typeStr, symbStr) 

150 

151 # -------------------------------------------------------------------------- 

152 # Properties. 

153 # -------------------------------------------------------------------------- 

154 

155 @property 

156 def dim(self): 

157 """The dimensionality :math:`n`.""" 

158 return self._n 

159 

160 @property 

161 def A(self): 

162 """The linear operator matrix :math:`A`.""" 

163 return self._A 

164 

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 

172 

173 return Constant(glyphs.inverse(self._A.string), inverse, self._A.shape) 

174 

175 @property 

176 def c(self): 

177 """The center point :math:`c`.""" 

178 return self._c 

179 

180 # -------------------------------------------------------------------------- 

181 # Abstract method implementations. 

182 # -------------------------------------------------------------------------- 

183 

184 def _get_mutables(self): 

185 return frozenset() 

186 

187 def _replace_mutables(self, mapping): 

188 return self 

189 

190 Subtype = namedtuple("Subtype", ("dim",)) 

191 

192 def _get_subtype(self): 

193 return self.Subtype(self._n) 

194 

195 # -------------------------------------------------------------------------- 

196 # Algebraic operations. 

197 # -------------------------------------------------------------------------- 

198 

199 @convert_operands() 

200 def __add__(self, other): 

201 if not isinstance(other, AffineExpression): 

202 return NotImplemented 

203 

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

205 

206 @convert_operands() 

207 def __radd__(self, other): 

208 if not isinstance(other, AffineExpression): 

209 return NotImplemented 

210 

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

212 

213 @convert_operands() 

214 def __mul__(self, other): 

215 if not isinstance(other, AffineExpression): 

216 return NotImplemented 

217 

218 if not other.scalar: 

219 raise TypeError( 

220 "Can only multiply an Ellipsoid from the right with a scalar.") 

221 

222 return Ellipsoid(self._n, self._A*other, self._c*other) 

223 

224 @convert_operands() 

225 def __rmul__(self, other): 

226 if not isinstance(other, AffineExpression): 

227 return NotImplemented 

228 

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

233 

234 return Ellipsoid(self._n, other*self._A, other*self._c) 

235 

236 @convert_operands() 

237 def __truediv__(self, other): 

238 if not isinstance(other, AffineExpression): 

239 return NotImplemented 

240 

241 if not other.scalar: 

242 raise TypeError("You may only divide an Ellipsoid by a scalar.") 

243 

244 if not other.constant: 

245 raise TypeError("You may only divide an Ellipsoid by a constant.") 

246 

247 if other.is0: 

248 raise ZeroDivisionError("Tried to divide an Ellipsoid by zero.") 

249 

250 return Ellipsoid(self._n, self._A/other, self._c/other) 

251 

252 # -------------------------------------------------------------------------- 

253 # Constraint-creating operations. 

254 # -------------------------------------------------------------------------- 

255 

256 @classmethod 

257 def _predict(cls, subtype, relation, other): 

258 assert isinstance(subtype, cls.Subtype) 

259 

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) 

264 

265 return NotImplemented 

266 

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

273 

274 one = AffineExpression.from_constant(1) 

275 return SOCConstraint(self.Ainv*(element.vec - self.c), one) 

276 else: 

277 return NotImplemented 

278 

279 

280# -------------------------------------- 

281__all__ = api_end(_API_START, globals())