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) 2019-2020 Maximilian Stahlberg 

3# Based on the original picos.expressions module by Guillaume Sagnol. 

4# 

5# This file is part of PICOS. 

6# 

7# PICOS is free software: you can redistribute it and/or modify it under the 

8# terms of the GNU General Public License as published by the Free Software 

9# Foundation, either version 3 of the License, or (at your option) any later 

10# version. 

11# 

12# PICOS is distributed in the hope that it will be useful, but WITHOUT ANY 

13# WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR 

14# A PARTICULAR PURPOSE. See the GNU General Public License for more details. 

15# 

16# You should have received a copy of the GNU General Public License along with 

17# this program. If not, see <http://www.gnu.org/licenses/>. 

18# ------------------------------------------------------------------------------ 

19 

20"""Implements affine expression types.""" 

21 

22import operator 

23from collections import namedtuple 

24 

25import numpy 

26 

27from .. import glyphs 

28from ..apidoc import api_end, api_start 

29from ..caching import cached_property, cached_unary_operator 

30from ..constraints import (AffineConstraint, ComplexAffineConstraint, 

31 ComplexLMIConstraint, Constraint, LMIConstraint) 

32from .data import convert_operands, cvx2np, load_data, sparse_quadruple 

33from .exp_biaffine import BiaffineExpression 

34from .expression import (Expression, ExpressionType, refine_operands, 

35 validate_prediction) 

36 

37_API_START = api_start(globals()) 

38# ------------------------------- 

39 

40 

41class ComplexAffineExpression(BiaffineExpression): 

42 """A multidimensional (complex) affine expression. 

43 

44 Base class for the real :class:`AffineExpression`. 

45 """ 

46 

47 # -------------------------------------------------------------------------- 

48 # Abstract method implementations for Expression. 

49 # -------------------------------------------------------------------------- 

50 

51 Subtype = namedtuple("Subtype", ("shape", "constant", "nonneg")) 

52 Subtype.dim = property(lambda self: self.shape[0] * self.shape[1]) 

53 

54 def _get_subtype(self): 

55 """Implement :meth:`~.expression.Expression._get_subtype`.""" 

56 nonneg = self.constant and self.isreal \ 

57 and all(x >= 0 for x in self.value_as_matrix) 

58 

59 return self.Subtype(self._shape, self.constant, nonneg) 

60 

61 # -------------------------------------------------------------------------- 

62 # Method overridings for Expression. 

63 # -------------------------------------------------------------------------- 

64 

65 @cached_unary_operator 

66 def _get_refined(self): 

67 """Implement :meth:`~.expression.Expression._get_refined`.""" 

68 if self.isreal: 

69 return AffineExpression(self._symbStr, self._shape, self._coefs) 

70 else: 

71 return self 

72 

73 @convert_operands(sameShape=True, allowNone=True) 

74 def _set_value(self, value): 

75 """Override :meth:`~.expression.Expression._set_value`.""" 

76 if value is None: 

77 for var in self._linear_coefs: 

78 var.value = None 

79 return 

80 

81 # Since all variables are real-valued, prevent NumPy from finding 

82 # complex-valued solutions that do not actually work. 

83 (self.real // self.imag).renamed(self.string).value \ 

84 = (value.real // value.imag) 

85 

86 # -------------------------------------------------------------------------- 

87 # Abstract method implementations for BiaffineExpression. 

88 # -------------------------------------------------------------------------- 

89 

90 @classmethod 

91 def _get_bilinear_terms_allowed(cls): 

92 """Implement for :class:`~.exp_biaffine.BiaffineExpression`.""" 

93 return False 

94 

95 @classmethod 

96 def _get_parameters_allowed(cls): 

97 """Implement for :class:`~.exp_biaffine.BiaffineExpression`.""" 

98 return False 

99 

100 @classmethod 

101 def _get_basetype(cls): 

102 """Implement :meth:`~.exp_biaffine.BiaffineExpression._get_basetype`.""" 

103 return ComplexAffineExpression 

104 

105 @classmethod 

106 def _get_typecode(cls): 

107 """Implement :meth:`~.exp_biaffine.BiaffineExpression._get_typecode`.""" 

108 return "z" 

109 

110 # -------------------------------------------------------------------------- 

111 # Method overridings for BiaffineExpression: Binary operators. 

112 # -------------------------------------------------------------------------- 

113 

114 @convert_operands(sameShape=True) 

115 @refine_operands(stop_at_affine=True) 

116 def __or__(self, other): 

117 """Extend :meth:`~.exp_biaffine.BiaffineExpression.__or__`. 

118 

119 The scalar product of two nonconstant (complex) affine expressions is in 

120 general a :class:`~.exp_quadratic.QuadraticExpression`. 

121 """ 

122 from .exp_quadratic import QuadraticExpression 

123 from .exp_sqnorm import SquaredNorm 

124 

125 if isinstance(other, ComplexAffineExpression) \ 

126 and not self.constant and not other.constant: 

127 # Create a squared norm if possible. 

128 # NOTE: Must not check self.equals(other) here; see SquaredNorm. 

129 # TODO: Consider creating a helper function for __or__ that always 

130 # returns a QuadraticExpression instead of a SquaredNorm to be 

131 # used within SquaredNorm. Then equals would be possible here. 

132 if self is other: 

133 return SquaredNorm(self) 

134 

135 string = glyphs.clever_dotp( 

136 self.string, other.string, other.complex, self.scalar) 

137 

138 # Handle the complex case: Conjugate the right hand side. 

139 other = other.conj 

140 

141 Cs, Co = self._constant_coef, other._constant_coef 

142 

143 # Compute the affine part of the product. 

144 affString = glyphs.affpart(string) 

145 affCoefs = {(): Cs.T * Co} 

146 for var in self.variables.union(other.variables): 

147 if var not in other._linear_coefs: 

148 affCoefs[var] = Co.T * self._linear_coefs[var] 

149 elif var not in self._linear_coefs: 

150 affCoefs[var] = Cs.T * other._linear_coefs[var] 

151 else: 

152 affCoefs[var] = Co.T * self._linear_coefs[var] + \ 

153 Cs.T * other._linear_coefs[var] 

154 affPart = self._common_basetype(other)(affString, (1, 1), affCoefs) 

155 

156 # Compute the quadratic part of the product. 

157 quadPart = {(v, w): self._linear_coefs[v].T * other._linear_coefs[w] 

158 for v in self._linear_coefs for w in other._linear_coefs} 

159 

160 # Don't create quadratic expressions without a quadratic part. 

161 if not any(quadPart.values()): 

162 affPart._symbStr = string 

163 return affPart 

164 

165 # Remember a factorization into two real scalars if applicable. 

166 # NOTE: If the user enters a multiplication a*b of two scalar affine 

167 # expressions, then we have, at this point, self == a.T == a 

168 # and other == b.conj.conj == b. 

169 if len(self) == 1 and len(other) == 1 \ 

170 and self.isreal and other.isreal: 

171 factors = (self.refined, other.refined) 

172 else: 

173 factors = None 

174 

175 return QuadraticExpression( 

176 string, quadPart, affPart, scalarFactors=factors) 

177 else: 

178 return BiaffineExpression.__or__(self, other) 

179 

180 @convert_operands(rMatMul=True) 

181 @refine_operands(stop_at_affine=True) 

182 def __mul__(self, other): 

183 """Extend :meth:`~.exp_biaffine.BiaffineExpression.__mul__`. 

184 

185 The matrix product of two nonconstant scalar (complex) affine 

186 expressions is a :class:`~.exp_quadratic.QuadraticExpression`. 

187 """ 

188 if isinstance(other, ComplexAffineExpression) \ 

189 and not self.constant and not other.constant: 

190 # If the result is scalar, allow for quadratic terms. 

191 if self._shape[0] == 1 and other._shape[1] == 1 \ 

192 and self._shape[1] == other._shape[0]: 

193 result = self.T.__or__(other.conj) 

194 

195 # NOTE: __or__ always creates a fresh expression. 

196 result._symbStr = glyphs.clever_mul(self.string, other.string) 

197 

198 return result 

199 else: 

200 raise NotImplementedError( 

201 "PICOS does not support multidimensional quadratic " 

202 "expressions at this point. More precisely, one factor must" 

203 " be constant or the result must be scalar.") 

204 else: 

205 return BiaffineExpression.__mul__(self, other) 

206 

207 @convert_operands(sameShape=True) 

208 @refine_operands(stop_at_affine=True) 

209 def __xor__(self, other): 

210 """Extend :meth:`~.exp_biaffine.BiaffineExpression.__xor__`. 

211 

212 The hadamard product of two nonconstant scalar (complex) affine 

213 expressions is a :class:`~.exp_quadratic.QuadraticExpression`. 

214 """ 

215 if isinstance(other, ComplexAffineExpression) \ 

216 and not self.constant and not other.constant: 

217 # If the result is scalar, allow for quadratic terms. 

218 if self._shape == (1, 1): 

219 result = self.__or__(other.conj) 

220 

221 # NOTE: __or__ always creates a fresh expression. 

222 result._symbStr = glyphs.hadamard(self.string, other.string) 

223 

224 return result 

225 else: 

226 raise NotImplementedError( 

227 "PICOS does not support multidimensional quadratic " 

228 "expressions at this point. More precisely, one factor must" 

229 " be constant or the result must be scalar.") 

230 else: 

231 return BiaffineExpression.__xor__(self, other) 

232 

233 # TODO: Create a quadratic expression from a scalar Kronecker prod. 

234 

235 # -------------------------------------------------------------------------- 

236 # Method overridings for BiaffineExpression: Unary operators. 

237 # -------------------------------------------------------------------------- 

238 

239 @cached_property 

240 def real(self): 

241 """Override :meth:`~.exp_biaffine.BiaffineExpression.real`. 

242 

243 The result is returned as an :meth:`AffineExpression`. 

244 """ 

245 return AffineExpression(glyphs.real(self.string), self._shape, 

246 {vars: coef.real() for vars, coef in self._coefs.items()}) 

247 

248 @cached_property 

249 def imag(self): 

250 """Override :meth:`~.exp_biaffine.BiaffineExpression.imag`. 

251 

252 The result is returned as an :meth:`AffineExpression`. 

253 """ 

254 return AffineExpression(glyphs.imag(self.string), self._shape, 

255 {vars: coef.imag() for vars, coef in self._coefs.items()}) 

256 

257 # -------------------------------------------------------------------------- 

258 # Additional unary operators. 

259 # -------------------------------------------------------------------------- 

260 

261 @cached_unary_operator 

262 def __abs__(self): 

263 """Denote the Euclidean or Frobenius norm of the expression.""" 

264 from . import Norm 

265 

266 return Norm(self) 

267 

268 # -------------------------------------------------------------------------- 

269 # Constraint-creating operators, and _predict. 

270 # -------------------------------------------------------------------------- 

271 

272 @classmethod 

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

274 assert isinstance(subtype, cls.Subtype) 

275 

276 from .set import Set 

277 

278 if relation == operator.__eq__: 

279 if issubclass(other.clstype, ComplexAffineExpression): 

280 return ComplexAffineConstraint.make_type(dim=subtype.dim) 

281 elif relation == operator.__lshift__: 

282 if issubclass(other.clstype, ComplexAffineExpression): 

283 return ComplexLMIConstraint.make_type(int(subtype.dim**0.5)) 

284 elif issubclass(other.clstype, Set): 

285 other_type = ExpressionType(cls, subtype) 

286 return other.predict(operator.__rshift__, other_type) 

287 elif relation == operator.__rshift__: 

288 if issubclass(other.clstype, ComplexAffineExpression): 

289 return ComplexLMIConstraint.make_type(int(subtype.dim**0.5)) 

290 

291 return NotImplemented 

292 

293 @convert_operands(sameShape=True) 

294 @validate_prediction 

295 @refine_operands() 

296 def __eq__(self, other): 

297 if isinstance(other, ComplexAffineExpression): 

298 return ComplexAffineConstraint(self, other) 

299 else: 

300 return NotImplemented 

301 

302 # Since we define __eq__, __hash__ is not inherited. Do this manually. 

303 __hash__ = Expression.__hash__ 

304 

305 def _lshift_implementation(self, other): 

306 if isinstance(other, ComplexAffineExpression): 

307 return ComplexLMIConstraint(self, Constraint.LE, other) 

308 else: 

309 return NotImplemented 

310 

311 def _rshift_implementation(self, other): 

312 if isinstance(other, ComplexAffineExpression): 

313 return ComplexLMIConstraint(self, Constraint.GE, other) 

314 else: 

315 return NotImplemented 

316 

317 # -------------------------------------------------------------------------- 

318 # Interface for PICOS-internal use. 

319 # -------------------------------------------------------------------------- 

320 

321 def sparse_rows( 

322 self, varOffsetMap, lowerTriangle=False, upperTriangle=False, 

323 indexFunction=None): 

324 """Return a sparse list representation of the expression. 

325 

326 The method is intended for internal use: It simplifies passing affine 

327 constraints to solvers that support only scalar constraints. The idea is 

328 to pose the constraint as a single (multidimensional) affine expression 

329 bounded by zero, and use the coefficients and the constant term of this 

330 expression to fill the solver's constraint matrix (with columns 

331 representing scalar variables and rows representing scalar constraints). 

332 

333 :param dict varOffsetMap: Maps variables to column offsets. 

334 :param bool lowerTriangle: Whether to return only the lower triangular 

335 part of the expression. 

336 :param bool upperTriangle: Whether to return only the upper triangular 

337 part of the expression. 

338 :param indexFunction: Instead of adding the local variable index to the 

339 value returned by varOffsetMap, use the return value of this 

340 function, that takes as argument the variable and its local index, 

341 as the "column index", which need not be an integer. When this 

342 parameter is passed, the parameter varOffsetMap is ignored. 

343 

344 :returns: A list of triples (J, V, c) where J contains column indices 

345 (representing scalar variables), V contains coefficients for each 

346 column index and c is a constant term. 

347 """ 

348 if lowerTriangle and upperTriangle: 

349 lowerTriangle = False 

350 upperTriangle = False 

351 

352 rows = [] 

353 numRows = len(self) 

354 m = self.size[0] 

355 

356 for row in range(numRows): 

357 i, j = row % m, row // m 

358 

359 if lowerTriangle and i < j: 

360 rows.append(None) 

361 elif upperTriangle and j < i: 

362 rows.append(None) 

363 else: 

364 rows.append([[], [], 0.0]) 

365 

366 for var, coef in self._linear_coefs.items(): 

367 V, I, J, _ = sparse_quadruple(coef) 

368 

369 for localCoef, localConIndex, localVarIndex in zip(V, I, J): 

370 row = rows[localConIndex] 

371 

372 if not row: 

373 continue 

374 

375 # TODO: Use a single parameter for both types. 

376 if indexFunction: 

377 row[0].append(indexFunction(var, localVarIndex)) 

378 else: 

379 row[0].append(varOffsetMap[var] + localVarIndex) 

380 

381 row[1].append(localCoef) 

382 

383 for localConIndex in range(numRows): 

384 row = rows[localConIndex] 

385 

386 if not row: 

387 continue 

388 

389 row[2] = self._constant_coef[localConIndex] 

390 

391 return [row for row in rows if row] 

392 

393 

394class AffineExpression(ComplexAffineExpression): 

395 """A multidimensional real affine expression.""" 

396 

397 # -------------------------------------------------------------------------- 

398 # Method overridings for BiaffineExpression. 

399 # -------------------------------------------------------------------------- 

400 

401 @property 

402 def isreal(self): 

403 """Always true for :class:`AffineExpression` instances.""" # noqa 

404 return True 

405 

406 @property 

407 def real(self): 

408 """The :class:`AffineExpression` as is.""" # noqa 

409 return self 

410 

411 @cached_property 

412 def imag(self): 

413 """A zero of same shape as the :class:`AffineExpression`.""" # noqa 

414 return self._basetype.zero(self._shape) 

415 

416 @property 

417 def conj(self): 

418 """The :class:`AffineExpression` as is.""" # noqa 

419 return self 

420 

421 @property 

422 def H(self): 

423 """The regular transpose of the :class:`AffineExpression`.""" # noqa 

424 return self.T 

425 

426 # -------------------------------------------------------------------------- 

427 # Method overridings for ComplexAffineExpression. 

428 # -------------------------------------------------------------------------- 

429 

430 @classmethod 

431 def _get_basetype(cls): 

432 return AffineExpression 

433 

434 @classmethod 

435 def _get_typecode(cls): 

436 return "d" 

437 

438 def _get_refined(self): 

439 return self 

440 

441 @convert_operands(sameShape=True, allowNone=True) 

442 def _set_value(self, value): 

443 if value is None: 

444 for var in self._linear_coefs: 

445 var.value = None 

446 return 

447 

448 if not isinstance(value, AffineExpression) or not value.constant: 

449 raise TypeError("Cannot set the value of {} to {}: Not real or not " 

450 "a constant.".format(repr(self), repr(value))) 

451 

452 if self.constant: 

453 raise TypeError("Cannot set the value on a constant expression.") 

454 

455 y = cvx2np(value._constant_coef) 

456 

457 A = [] 

458 for var, coef in self._linear_coefs.items(): 

459 A.append(cvx2np(coef)) 

460 assert A 

461 

462 A = numpy.hstack(A) 

463 b = y - cvx2np(self._constant_coef) 

464 

465 try: 

466 solution, residual, _, _ = numpy.linalg.lstsq(A, b, rcond=None) 

467 except numpy.linalg.LinAlgError as error: 

468 raise RuntimeError("Setting a value on {} by means of a least-" 

469 "squares solution failed.".format(self.string)) from error 

470 

471 if not numpy.allclose(residual, 0): 

472 raise ValueError("Setting a value on {} failed: No exact solution " 

473 "to the associated linear system found.".format(self.string)) 

474 

475 offset = 0 

476 for var in self._linear_coefs: 

477 var.internal_value = solution[offset:offset+var.dim] 

478 offset += var.dim 

479 

480 # -------------------------------------------------------------------------- 

481 # Additional unary operators. 

482 # -------------------------------------------------------------------------- 

483 

484 @cached_property 

485 def exp(self): 

486 """The exponential function applied to the expression.""" # noqa 

487 from . import SumExponentials 

488 return SumExponentials(self) 

489 

490 @cached_property 

491 def log(self): 

492 """The Logarithm of the expression.""" # noqa 

493 from . import Logarithm 

494 return Logarithm(self) 

495 

496 # -------------------------------------------------------------------------- 

497 # Constraint-creating operators, and _predict. 

498 # -------------------------------------------------------------------------- 

499 

500 @classmethod 

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

502 assert isinstance(subtype, cls.Subtype) 

503 

504 if relation in (operator.__eq__, operator.__le__, operator.__ge__): 

505 if issubclass(other.clstype, AffineExpression): 

506 return AffineConstraint.make_type( 

507 dim=subtype.dim, eq=(relation is operator.__eq__)) 

508 elif relation == operator.__lshift__: 

509 if issubclass(other.clstype, AffineExpression): 

510 return LMIConstraint.make_type(int(subtype.dim**0.5)) 

511 elif issubclass(other.clstype, ComplexAffineExpression): 

512 return ComplexLMIConstraint.make_type(int(subtype.dim**0.5)) 

513 elif relation == operator.__rshift__: 

514 if issubclass(other.clstype, AffineExpression): 

515 return LMIConstraint.make_type(int(subtype.dim**0.5)) 

516 elif issubclass(other.clstype, ComplexAffineExpression): 

517 return ComplexLMIConstraint.make_type(int(subtype.dim**0.5)) 

518 

519 return NotImplemented 

520 

521 @convert_operands(sameShape=True) 

522 @validate_prediction 

523 @refine_operands() 

524 def __le__(self, other): 

525 if isinstance(other, AffineExpression): 

526 return AffineConstraint(self, Constraint.LE, other) 

527 else: 

528 return NotImplemented 

529 

530 @convert_operands(sameShape=True) 

531 @validate_prediction 

532 @refine_operands() 

533 def __ge__(self, other): 

534 if isinstance(other, AffineExpression): 

535 return AffineConstraint(self, Constraint.GE, other) 

536 else: 

537 return NotImplemented 

538 

539 @convert_operands(sameShape=True) 

540 @validate_prediction 

541 @refine_operands() 

542 def __eq__(self, other): 

543 if isinstance(other, AffineExpression): 

544 return AffineConstraint(self, Constraint.EQ, other) 

545 else: 

546 return NotImplemented 

547 

548 # Since we define __eq__, __hash__ is not inherited. Do this manually. 

549 __hash__ = Expression.__hash__ 

550 

551 def _lshift_implementation(self, other): 

552 if isinstance(other, AffineExpression): 

553 return LMIConstraint(self, Constraint.LE, other) 

554 elif isinstance(other, ComplexAffineExpression): 

555 return ComplexLMIConstraint(self, Constraint.LE, other) 

556 else: 

557 return NotImplemented 

558 

559 def _rshift_implementation(self, other): 

560 if isinstance(other, AffineExpression): 

561 return LMIConstraint(self, Constraint.GE, other) 

562 elif isinstance(other, ComplexAffineExpression): 

563 return ComplexLMIConstraint(self, Constraint.GE, other) 

564 else: 

565 return NotImplemented 

566 

567 

568def Constant(name_or_value, value=None, shape=None): 

569 """Create a constant PICOS expression. 

570 

571 Loads the given numeric value as a constant 

572 :class:`~picos.expressions.ComplexAffineExpression` or 

573 :class:`~picos.expressions.AffineExpression`, depending on the value. 

574 Optionally, the value is broadcasted or reshaped according to the shape 

575 argument. 

576 

577 :param str name_or_value: Symbolic string description of the constant. If 

578 :obj:`None` or the empty string, a string will be generated. If this is 

579 the only positional parameter (i.e.``value`` is not given), then this 

580 position is used as the value argument instead! 

581 :param value: The numeric constant to load. 

582 

583 See :func:`~.data.load_data` for supported data formats and broadcasting and 

584 reshaping rules. 

585 

586 :Example: 

587 

588 >>> from picos import Constant 

589 >>> Constant(1) 

590 <1×1 Real Constant: 1> 

591 >>> Constant(1, shape=(2, 2)) 

592 <2×2 Real Constant: [1]> 

593 >>> Constant("one", 1) 

594 <1×1 Real Constant: one> 

595 >>> Constant("J", 1, (2, 2)) 

596 <2×2 Real Constant: J> 

597 """ 

598 if value is None: 

599 value = name_or_value 

600 name = None 

601 else: 

602 name = name_or_value 

603 

604 value, valStr = load_data(value, shape) 

605 

606 if value.typecode == "z": 

607 cls = ComplexAffineExpression 

608 else: 

609 cls = AffineExpression 

610 

611 return cls(name if name else valStr, value.size, {(): value}) 

612 

613 

614# -------------------------------------- 

615__all__ = api_end(_API_START, globals())