Coverage for picos/expressions/set

1# ------------------------------------------------------------------------------

4# This file is part of PICOS.

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

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

Coverage for picos/expressions/set_ellipsoid.py: 69.91%

113 statements