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 :class:`BiaffineExpression`.""" 

21 

22import math 

23import operator 

24from abc import ABC, abstractmethod 

25from functools import reduce 

26from itertools import product 

27 

28import cvxopt 

29import numpy 

30 

31from .. import glyphs, settings 

32from ..apidoc import api_end, api_start 

33from ..caching import (cached_property, cached_selfinverse_property, 

34 cached_selfinverse_unary_operator, 

35 cached_unary_operator, unlocked_cached_properties) 

36from ..formatting import detect_range 

37from ..legacy import deprecated 

38from .data import (blend_shapes, convert_operands, cvx2np, cvxopt_equals, 

39 cvxopt_K, cvxopt_vcat, left_kronecker_I, load_data, 

40 load_dense_data, load_shape, right_kronecker_I) 

41from .expression import Expression, refine_operands 

42from .mutable import Mutable 

43from .vectorizations import FullVectorization, SymmetricVectorization 

44 

45_API_START = api_start(globals()) 

46# ------------------------------- 

47 

48 

49class BiaffineExpression(Expression, ABC): 

50 r"""A multidimensional (complex) biaffine expression. 

51 

52 Abstract base class for the affine 

53 :class:`~.exp_affine.ComplexAffineExpression` and 

54 its real subclass :class:`~.exp_affine.AffineExpression`. 

55 Quadratic expressions are stored in 

56 :class:`~.exp_quadratic.QuadraticExpression` instead. 

57 

58 In general this expression has the form 

59 

60 .. math:: 

61 

62 A(x,y) = B(x,y) + P(x) + Q(y) + C 

63 

64 where :math:`x \in \mathbb{R}^p` and :math:`y \in \mathbb{R}^q` are variable 

65 vectors, 

66 :math:`B : \mathbb{R}^p \times \mathbb{R}^q \to \mathbb{C}^{m \times n}` 

67 is a bilinear function, 

68 :math:`P : \mathbb{R}^p \to \mathbb{C}^{m \times n}` and 

69 :math:`Q : \mathbb{R}^q \to \mathbb{C}^{m \times n}` are linear functions, 

70 and :math:`C \in \mathbb{C}^{m \times n}` is a constant. 

71 

72 If no coefficient matrices defining :math:`B` and :math:`Q` are provided on 

73 subclass instanciation, then this acts as an affine function of :math:`x`. 

74 

75 In a more technical sense, the notational variables :math:`x` and :math:`y` 

76 each represent a stack of vectorizations of a number of *actual* scalar, 

77 vector, or matrix variables or parameters :math:`X_i` and :math:`Y_j` with 

78 :math:`i, j \in \mathbb{Z}_{\geq 0}` and this class stores the nonzero 

79 (bi)linear functions 

80 :math:`B_{i,j}(\operatorname{vec}(X_i),\operatorname{vec}(Y_j))`, 

81 :math:`P_{i}(\operatorname{vec}(X_i))` and 

82 :math:`Q_{j}(\operatorname{vec}(Y_j))` in the form of separate sparse 

83 coefficient matrices. 

84 """ 

85 

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

87 # Static methods. 

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

89 

90 @staticmethod 

91 def _update_coefs(coefs, mtbs, summand): 

92 if not summand: 

93 return 

94 elif mtbs in coefs: 

95 coefs[mtbs] = coefs[mtbs] + summand 

96 else: 

97 coefs[mtbs] = summand 

98 

99 # -------------------------------------------------------------------------- 

100 # Initialization and factory methods. 

101 # -------------------------------------------------------------------------- 

102 

103 @classmethod 

104 def _get_type_string_base(cls): 

105 """Return a string template for the expression type string.""" 

106 return "Complex {}" if cls._get_typecode() == "z" else "Real {}" 

107 

108 def __init__(self, string, shape=(1, 1), coefficients={}): 

109 """Initialize a (complex) biaffine expression. 

110 

111 :param str string: A symbolic string description. 

112 

113 :param shape: Shape of a vector or matrix expression. 

114 :type shape: int or tuple or list 

115 

116 :param dict coefficients: 

117 Maps mutable pairs, single mutables or the empty tuple to sparse 

118 represetations of the coefficient matrices :math:`B`, :math:`P` and 

119 :math:`Q`, and :math:`C`, respectively. 

120 

121 .. warning:: 

122 If the given coefficients are already of the desired (numeric) type 

123 and shape, they are stored by reference. Modifying such data can 

124 lead to unexpected results as PICOS expressions are supposed to be 

125 immutable (to allow 

126 `caching <https://en.wikipedia.org/wiki/Cache_(computing)>`_ of 

127 results). 

128 

129 If you create biaffine expressions by any other means than this 

130 constructor, PICOS makes a copy of your data to prevent future 

131 modifications to it from causing inconsistencies. 

132 """ 

133 from .variables import BaseVariable 

134 

135 shape = load_shape(shape) 

136 length = shape[0]*shape[1] 

137 

138 if not isinstance(coefficients, dict): 

139 raise TypeError("Coefficients of {} must be given as a dict." 

140 .format(type(self).__name__)) 

141 

142 # Store shape. 

143 self._shape = shape 

144 

145 # Store coefficients. 

146 self._coefs = {} 

147 for mtbs, coef in coefficients.items(): 

148 if isinstance(mtbs, Mutable): 

149 mtbs = (mtbs,) 

150 elif not isinstance(mtbs, tuple): 

151 raise TypeError("Coefficient indices must be tuples, not " 

152 "objects of type {}.".format(type(mtbs).__name__)) 

153 

154 if not all(isinstance(mtb, Mutable) for mtb in mtbs): 

155 raise TypeError("Coefficients must be indexed by mutables.") 

156 

157 if not self._parameters_allowed \ 

158 and not all(isinstance(mtb, BaseVariable) for mtb in mtbs): 

159 raise TypeError("Coefficients of {} may only be indexed by " 

160 "decision variables.".format(type(self).__name__)) 

161 

162 # Store only linear and potentially bilinear terms. 

163 if len(mtbs) != len(set(mtbs)): 

164 raise TypeError("Coefficients of {} must be indexed by disjoint" 

165 " mutables; no quadratic terms are allowed." 

166 .format(type(self).__name__)) 

167 elif len(mtbs) > 2: 

168 raise TypeError("Coefficients of {} may be indexed by at " 

169 "most two mutables.".format(type(self).__name__)) 

170 elif len(mtbs) > 1 and not self._bilinear_terms_allowed: 

171 raise TypeError("Coefficients of {} may be indexed by at " 

172 "most one mutable.".format(type(self).__name__)) 

173 

174 # Do not store coefficients of zero. 

175 if not coef: 

176 continue 

177 

178 # Dimension of the tensor product of the vectorized mutables. 

179 dim = reduce(operator.mul, (mtb.dim for mtb in mtbs), 1) 

180 

181 # Load the coefficient matrix in the desired format. 

182 coef = load_data( 

183 coef, (length, dim), self._typecode, alwaysCopy=False)[0] 

184 

185 # Always store with respect to ordered mutables. 

186 if len(mtbs) == 2 and mtbs[0].id > mtbs[1].id: 

187 # Obtain a fitting commutation matrix. 

188 K = cvxopt_K(mtbs[1].dim, mtbs[0].dim, self._typecode) 

189 

190 # Make coef apply to vec(y*x.T) instead of vec(x*y.T). 

191 coef = coef * K 

192 

193 # Swap x and y. 

194 mtbs = (mtbs[1], mtbs[0]) 

195 

196 # Sum with existing coefficients submitted in opposing order. 

197 if mtbs in self._coefs: 

198 self._coefs[mtbs] = self._coefs[mtbs] + coef 

199 else: 

200 self._coefs[mtbs] = coef 

201 

202 # Determine the type string. 

203 typeStr = self._get_type_string_base() 

204 

205 if "{}" in typeStr: 

206 hasBilinear = bool(self._bilinear_coefs) 

207 hasLinear = bool(self._linear_coefs) 

208 hasConstant = bool(self._constant_coef) 

209 

210 if hasBilinear: 

211 if hasConstant: 

212 typeStr = typeStr.format("Biaffine Expression") 

213 else: 

214 typeStr = typeStr.format("Bilinear Expression") 

215 elif hasLinear: 

216 if hasConstant: 

217 typeStr = typeStr.format("Affine Expression") 

218 else: 

219 typeStr = typeStr.format("Linear Expression") 

220 else: 

221 typeStr = typeStr.format("Constant") 

222 

223 typeStr = "{} {}".format(glyphs.shape(shape), typeStr) 

224 

225 Expression.__init__(self, typeStr, string) 

226 

227 @classmethod 

228 def from_constant(cls, constant, shape=None, name=None): 

229 """Create a class instance from the given numeric constant. 

230 

231 Loads the given constant as a PICOS expression, optionally broadcasted 

232 or reshaped to the given shape and named as specified. 

233 

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

235 and reshaping rules. 

236 

237 Unlike :func:`~.exp_affine.Constant`, this class method always creates 

238 an instance of the class that it is called on, instead of tailoring 

239 towards the numeric type of the data. 

240 

241 .. note:: 

242 When an operation involves both a PICOS expression and a constant 

243 value of another type, PICOS converts the constant on the fly so 

244 that you rarely need to use this method. 

245 """ 

246 constant, string = load_data(constant, shape, cls._get_typecode()) 

247 return cls._get_basetype()( 

248 name if name else string, constant.size, {(): constant[:]}) 

249 

250 @classmethod 

251 def zero(cls, shape=(1, 1)): 

252 """Return a constant zero expression of given shape.""" 

253 shape = load_shape(shape) 

254 string = glyphs.scalar(0) if shape == (1, 1) else glyphs.matrix(0) 

255 return cls._get_basetype()(string, shape) 

256 

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

258 # (Abstract) class methods. 

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

260 

261 @classmethod 

262 @abstractmethod 

263 def _get_bilinear_terms_allowed(cls): 

264 """Report whether the subclass may have bilinear terms.""" 

265 pass 

266 

267 @classmethod 

268 @abstractmethod 

269 def _get_parameters_allowed(cls): 

270 """Report whether the subclass may depend on parameters.""" 

271 pass 

272 

273 @classmethod 

274 @abstractmethod 

275 def _get_basetype(cls): 

276 """Return first non-abstract :class:`Expression` subclass this bases on. 

277 

278 Enables subclass objects (such as variables) to behave like the returned 

279 type with respect to algebraic operations. For instance, the sum of 

280 two :class:`ComplexVariable` is a :class:`ComplexAffineExpression`. 

281 """ 

282 pass 

283 

284 @classmethod 

285 @abstractmethod 

286 def _get_typecode(cls): 

287 """Return the CVXOPT typecode to use with coefficient matrices. 

288 

289 Either ``"z"`` for complex or ``"d"`` for real. 

290 

291 See also :meth:`_get_basetype`. 

292 """ 

293 pass 

294 

295 @classmethod 

296 def _common_basetype(cls, other, reverse=False): 

297 """Return the common basetype of two expressions. 

298 

299 The parameter ``other`` may be either a :class:`BiaffineExpression` 

300 instance or subclass thereof. 

301 """ 

302 myBasetype = cls._get_basetype() 

303 theirBasetype = other._get_basetype() 

304 

305 if myBasetype is theirBasetype: 

306 return myBasetype 

307 elif issubclass(theirBasetype, myBasetype): 

308 return myBasetype 

309 elif issubclass(myBasetype, theirBasetype): 

310 return theirBasetype 

311 elif not reverse: 

312 # HACK: Handle the case where the individual base types do not have 

313 # a sub- and superclass relation but where one of the types 

314 # still knows what the resulting base class should be. 

315 return other._common_basetype(cls, reverse=True) 

316 else: 

317 raise TypeError("The expression types {} and {} do not have a " 

318 "common base type apart from the abstract {} so the operation " 

319 "cannot be performed.".format(cls.__name__, other.__name__, 

320 BiaffineExpression.__name__)) 

321 

322 # -------------------------------------------------------------------------- 

323 # Internal use properties and methods. 

324 # -------------------------------------------------------------------------- 

325 

326 @property 

327 def _bilinear_terms_allowed(self): 

328 """Shorthand for :meth:`_get_bilinear_terms_allowed`.""" 

329 return self._get_bilinear_terms_allowed() 

330 

331 @property 

332 def _parameters_allowed(self): 

333 """Shorthand for :meth:`_get_parameters_allowed`.""" 

334 return self._get_parameters_allowed() 

335 

336 @property 

337 def _basetype(self): 

338 """Shorthand for :meth:`_get_basetype`.""" 

339 return self._get_basetype() 

340 

341 @property 

342 def _typecode(self): 

343 """Shorthand for :meth:`_get_typecode`.""" 

344 return self._get_typecode() 

345 

346 @cached_property 

347 def _constant_coef(self): 

348 """Vectorized constant term with zero represented explicitly.""" 

349 if () in self._coefs: 

350 return self._coefs[()] 

351 else: 

352 return load_data(0, len(self), self._typecode)[0] 

353 

354 @cached_property 

355 def _linear_coefs(self): 

356 """Linear part coefficients indexed directly by single mutables.""" 

357 return {mtbs[0]: c for mtbs, c in self._coefs.items() if len(mtbs) == 1} 

358 

359 @cached_property 

360 def _bilinear_coefs(self): 

361 """Bilinear part coefficients indexed by mutable pairs.""" 

362 return {mtbs: c for mtbs, c in self._coefs.items() if len(mtbs) == 2} 

363 

364 @cached_property 

365 def _sparse_coefs(self): 

366 """Coefficients (re-)cast to sparse matrices.""" 

367 return {mtbs: cvxopt.sparse(coef) for mtbs, coef in self._coefs.items()} 

368 

369 @cached_property 

370 def _sparse_linear_coefs(self): 

371 """Linear part coefficients cast to sparse and indexed by mutables.""" 

372 return {v[0]: c for v, c in self._sparse_coefs.items() if len(v) == 1} 

373 

374 def _reglyphed_string(self, mutable_name_map): 

375 """The symbolic string with mutable names replaced.""" 

376 string = self.string 

377 

378 if isinstance(string, glyphs.GlStr): 

379 string = string.reglyphed(mutable_name_map) 

380 elif string in mutable_name_map: 

381 # Handle corner cases like x + 0 for a mutable x which is not a 

382 # mutable any more but still has a literal string name. 

383 string = mutable_name_map[string] 

384 

385 return string 

386 

387 # -------------------------------------------------------------------------- 

388 # Abstract method implementations and overridings for Expression. 

389 # -------------------------------------------------------------------------- 

390 

391 def _get_value(self): 

392 # Create a copy of the constant term. 

393 value = self._constant_coef[:] 

394 

395 # Add value of the linear part. 

396 for mtb, coef in self._linear_coefs.items(): 

397 summand = coef * mtb._get_internal_value() 

398 

399 if type(value) == type(summand): 

400 value += summand 

401 else: 

402 # Exactly one of the matrices is sparse. 

403 value = value + summand 

404 

405 # Add value of the bilinear part. 

406 for (x, y), coef in self._bilinear_coefs.items(): 

407 xValue = x._get_internal_value() 

408 yValue = y._get_internal_value() 

409 summand = coef * (xValue*yValue.T)[:] 

410 

411 if type(value) == type(summand): 

412 value += summand 

413 else: 

414 # Exactly one of the matrices is sparse. 

415 value = value + summand 

416 

417 # Resize the value to the proper shape. 

418 value.size = self.shape 

419 

420 return value 

421 

422 def _get_shape(self): 

423 return self._shape 

424 

425 @cached_unary_operator 

426 def _get_mutables(self): 

427 return frozenset(mtb for mtbs in self._coefs for mtb in mtbs) 

428 

429 def _is_convex(self): 

430 return not self._bilinear_coefs 

431 

432 def _is_concave(self): 

433 return not self._bilinear_coefs 

434 

435 def _replace_mutables(self, mapping): 

436 # Handle the base case where the affine expression is a mutable. 

437 if self in mapping: 

438 return mapping[self] 

439 

440 name_map = {old.name: new.string for old, new in mapping.items()} 

441 string = self._reglyphed_string(name_map) 

442 

443 if all(isinstance(new, Mutable) for new in mapping.values()): 

444 # Fast implementation for the basic case. 

445 coefs = {tuple(mapping[mtb] for mtb in mtbs): coef 

446 for mtbs, coef in self._coefs.items()} 

447 else: 

448 # Turn full mapping into an effective mapping. 

449 mapping = {o: n for o, n in mapping.items() if o is not n} 

450 

451 coefs = {} 

452 for old_mtbs, old_coef in self._coefs.items(): 

453 if not any(old_mtb in mapping for old_mtb in old_mtbs): 

454 self._update_coefs(coefs, old_mtbs, old_coef) 

455 elif len(old_mtbs) == 1: 

456 assert old_mtbs[0] in mapping 

457 new = mapping[old_mtbs[0]] 

458 for new_mtb, new_coef in new._linear_coefs.items(): 

459 self._update_coefs(coefs, (new_mtb,), old_coef*new_coef) 

460 self._update_coefs(coefs, (), old_coef*new._constant_coef) 

461 elif old_mtbs[0] in mapping and old_mtbs[1] not in mapping: 

462 new = mapping[old_mtbs[0]] 

463 old_mtb = old_mtbs[1] 

464 for new_mtb, new_coef in new._linear_coefs.items(): 

465 self._update_coefs(coefs, (new_mtb, old_mtb), 

466 old_coef*(left_kronecker_I(new_coef, old_mtb.dim))) 

467 self._update_coefs(coefs, (old_mtb,), old_coef*( 

468 left_kronecker_I(new._constant_coef, old_mtb.dim))) 

469 elif old_mtbs[0] not in mapping and old_mtbs[1] in mapping: 

470 new = mapping[old_mtbs[1]] 

471 old_mtb = old_mtbs[0] 

472 for new_mtb, new_coef in new._linear_coefs.items(): 

473 self._update_coefs(coefs, (old_mtb, new_mtb), 

474 old_coef*(right_kronecker_I(new_coef, old_mtb.dim))) 

475 self._update_coefs(coefs, (old_mtb,), old_coef*( 

476 right_kronecker_I(new._constant_coef, old_mtb.dim))) 

477 elif old_mtbs[0] in mapping and old_mtbs[1] in mapping: 

478 new1 = mapping[old_mtbs[0]] 

479 new2 = mapping[old_mtbs[1]] 

480 if isinstance(new1, Mutable) and isinstance(new2, Mutable): 

481 self._update_coefs(coefs, (new1, new2), old_coef) 

482 else: 

483 raise NotImplementedError( 

484 "Replacing both mutables in a bilinear term is not " 

485 "supported unless both are replaced with mutables. " 

486 "The effective mapping was: {}".format(mapping)) 

487 else: 

488 assert False 

489 

490 return self._basetype(string, self._shape, coefs) 

491 

492 def _freeze_mutables(self, freeze): 

493 string = self._reglyphed_string( 

494 {mtb.name: glyphs.frozen(mtb.name) for mtb in freeze}) 

495 

496 coefs = {} 

497 for mtbs, coef in self._coefs.items(): 

498 if not any(mtb in freeze for mtb in mtbs): 

499 self._update_coefs(coefs, mtbs, coef) 

500 elif len(mtbs) == 1: 

501 assert mtbs[0] in freeze 

502 self._update_coefs(coefs, (), coef*mtbs[0].internal_value) 

503 elif mtbs[0] in freeze and mtbs[1] in freeze: 

504 C = coef*(mtbs[0].internal_value*mtbs[1].internal_value.T)[:] 

505 self._update_coefs(coefs, (), C) 

506 elif mtbs[0] in freeze and mtbs[1] not in freeze: 

507 C = coef*left_kronecker_I(mtbs[0].internal_value, mtbs[1].dim) 

508 self._update_coefs(coefs, (mtbs[1],), C) 

509 elif mtbs[0] not in freeze and mtbs[1] in freeze: 

510 C = coef*right_kronecker_I(mtbs[1].internal_value, mtbs[0].dim) 

511 self._update_coefs(coefs, (mtbs[0],), C) 

512 else: 

513 assert False 

514 

515 return self._basetype(string, self._shape, coefs) 

516 

517 # -------------------------------------------------------------------------- 

518 # Python special method implementations. 

519 # -------------------------------------------------------------------------- 

520 

521 def __len__(self): 

522 # Faster version that overrides Expression.__len__. 

523 return self._shape[0] * self._shape[1] 

524 

525 def __getitem__(self, index): 

526 def slice2range(s, length): 

527 """Transform a :class:`slice` to a :class:`range`.""" 

528 assert isinstance(s, slice) 

529 

530 # Plug in slice's default values. 

531 ss = s.step if s.step else 1 

532 if ss > 0: 

533 sa = s.start if s.start is not None else 0 

534 sb = s.stop if s.stop is not None else length 

535 else: 

536 assert ss < 0 

537 sa = s.start if s.start is not None else length - 1 

538 sb = s.stop # Keep it None as -1 would mean length - 1. 

539 

540 # Wrap negative indices (once). 

541 ra = length + sa if sa < 0 else sa 

542 if sb is None: 

543 # This is the only case where we give a negative index to range. 

544 rb = -1 

545 else: 

546 rb = length + sb if sb < 0 else sb 

547 

548 # Clamp out-of-bound indices. 

549 ra = min(max(0, ra), length - 1) 

550 rb = min(max(-1, rb), length) 

551 

552 r = range(ra, rb, ss) 

553 

554 if not r: 

555 raise IndexError("Empty slice.") 

556 

557 return r 

558 

559 def range2slice(r, length): 

560 """Transform a :class:`range` to a :class:`slice`, if possible. 

561 

562 :raises ValueError: If the input cannot be expressed as a slice. 

563 """ 

564 assert isinstance(r, range) 

565 

566 if not r: 

567 raise IndexError("Empty range.") 

568 

569 ra = r.start 

570 rb = r.stop 

571 rs = r.step 

572 

573 if rs > 0: 

574 if ra < 0 or rb > length: 

575 raise ValueError( 

576 "Out-of-bounds range cannot be represented as a slice.") 

577 else: 

578 assert rs < 0 

579 if ra >= length or rb < -1: 

580 raise ValueError( 

581 "Out-of-bounds range cannot be represented as a slice.") 

582 

583 if rb == -1: 

584 rb = None 

585 

586 return slice(ra, rb, rs) 

587 

588 def list2slice(l, length): 

589 """Transform a :class:`list` to a :class:`slice`, if possible. 

590 

591 :raises TypeError: If the input is not an integer sequence. 

592 :raises ValueError: If the input cannot be expressed as a slice. 

593 """ 

594 return range2slice(detect_range(l), length) 

595 

596 def slice2str(s): 

597 """Return the short string that produced a :class:`slice`.""" 

598 assert isinstance(s, slice) 

599 

600 startStr = str(s.start) if s.start is not None else "" 

601 stopStr = str(s.stop) if s.stop is not None else "" 

602 if s.step in (None, 1): 

603 if s.start is not None and s.stop is not None \ 

604 and s.stop == s.start + 1: 

605 return startStr 

606 else: 

607 return "{}:{}".format(startStr, stopStr) 

608 else: 

609 return "{}:{}:{}".format(startStr, stopStr, str(s.step)) 

610 

611 def list2str(l, length): 

612 """Return a short string represnetation of a :class:`list`.""" 

613 assert isinstance(l, list) 

614 

615 # Extract integers wrapped in a list. 

616 if len(l) == 1: 

617 return str(l[0]) 

618 

619 # Try to convert the list to a slice. 

620 try: 

621 l = list2slice(l, length) 

622 except (ValueError, RuntimeError): 

623 pass 

624 

625 if isinstance(l, list): 

626 if len(l) > 4: 

627 return glyphs.shortint(l[0], l[-1]) 

628 else: 

629 return str(l).replace(" ", "") 

630 else: 

631 return slice2str(l) 

632 

633 def any2str(a, length): 

634 if isinstance(a, slice): 

635 return slice2str(a) 

636 elif isinstance(a, list): 

637 return list2str(a, length) 

638 else: 

639 assert False 

640 

641 m, n = self._shape 

642 indexStr = None 

643 isIntList = False 

644 

645 # Turn the index expression into a mutable list of one index per axis. 

646 if isinstance(index, tuple): # Multiple axis slicing. 

647 index = list(index) 

648 elif isinstance(index, dict): # Arbitrary matrix element selection. 

649 if len(index) != 2: 

650 raise TypeError("When slicing with a dictionary, there must be " 

651 "exactly two keys.") 

652 

653 try: 

654 i, j = sorted(index.keys()) 

655 except TypeError as error: 

656 raise TypeError("When slicing with a dictionary, the two keys " 

657 "must be comparable.") from error 

658 

659 I, J = index[i], index[j] 

660 

661 try: 

662 I = load_dense_data(I, typecode="i", alwaysCopy=False)[0] 

663 J = load_dense_data(J, typecode="i", alwaysCopy=False)[0] 

664 

665 if 1 not in I.size or 1 not in J.size: 

666 raise TypeError("At least one of the objects is not flat " 

667 "but a proper matrix.") 

668 

669 if len(I) != len(J): 

670 raise TypeError("The objects do not have the same length.") 

671 

672 I, J = list(I), list(J) 

673 except Exception as error: 

674 raise TypeError("Loading a sparse index vector pair for {} from" 

675 " objects of type {} and {} failed: {}".format(self.string, 

676 type(index[i]).__name__, type(index[j]).__name__, error)) \ 

677 from None 

678 

679 # Represent the selection as a global index list. 

680 index = [[i + j*m for i, j in zip(I, J)]] 

681 

682 # Use a special index string. 

683 indexStr = glyphs.size(list2str(I, n), list2str(J, m)) 

684 

685 # Don't invoke load_dense_data on the list. 

686 isIntList = True 

687 else: # Global indexing. 

688 index = [index] 

689 

690 # Make sure either global or row/column indexing is used. 

691 if not index: 

692 raise IndexError("Empty index.") 

693 elif len(index) > 2: 

694 raise IndexError( 

695 "PICOS expressions do not have a third axis to slice.") 

696 

697 # Turn the indices for each axis into either a slice or a list. 

698 for axis, selection in enumerate(index): 

699 # Convert anything that is not a slice, including scalars and lists 

700 # that are not confirmed integer, to an integer row or column 

701 # vector, then (back) to a list. 

702 if not isIntList and not isinstance(selection, slice): 

703 try: 

704 matrix = load_dense_data( 

705 selection, typecode="i", alwaysCopy=False)[0] 

706 

707 if 1 not in matrix.size: 

708 raise TypeError("The object is not flat but a {} shaped" 

709 " matrix.".format(glyphs.shape(matrix.size))) 

710 

711 selection = list(matrix) 

712 except Exception as error: 

713 raise TypeError("Loading a slicing index vector for axis {}" 

714 " of {} from an object of type {} failed: {}".format( 

715 axis, self.string, type(selection).__name__, error)) \ 

716 from None 

717 

718 index[axis] = selection 

719 

720 # Build index string, retrieve new shape, finalize index. 

721 if len(index) == 1: # Global indexing. 

722 index = index[0] 

723 

724 if isinstance(index, slice): 

725 shape = len(slice2range(index, len(self))) 

726 else: 

727 shape = len(index) 

728 

729 if indexStr is None: 

730 indexStr = any2str(index, len(self)) 

731 else: # Multiple axis slicing. 

732 if indexStr is None: 

733 indexStr = "{},{}".format( 

734 any2str(index[0], m), any2str(index[1], n)) 

735 

736 # Convert to a global index list. 

737 RC, shape = [], [] 

738 for axis, selection in enumerate(index): 

739 k = self._shape[axis] 

740 

741 if isinstance(selection, slice): 

742 # Turn the slice into an iterable range. 

743 selection = slice2range(selection, k) 

744 

745 # All indices from a slice are nonnegative. 

746 assert all(i >= 0 for i in selection) 

747 

748 if isinstance(selection, list): 

749 # Wrap once. This is consistent with CVXOPT. 

750 selection = [i if i >= 0 else k + i for i in selection] 

751 

752 # Perform a partial out-of-bounds check. 

753 if any(i < 0 for i in selection): 

754 raise IndexError( 

755 "Out-of-bounds access along axis {}.".format(axis)) 

756 

757 # Complete the check for out-of-bounds access. 

758 if any(i >= k for i in selection): 

759 raise IndexError( 

760 "Out-of-bounds access along axis {}.".format(axis)) 

761 

762 RC.append(selection) 

763 shape.append(len(selection)) 

764 

765 rows, cols = RC 

766 index = [i + j*m for j in cols for i in rows] 

767 

768 # Finalize the string. 

769 string = glyphs.slice(self.string, indexStr) 

770 

771 # Retrieve new coefficients and constant term. 

772 coefs = {mtbs: coef[index, :] for mtbs, coef in self._coefs.items()} 

773 

774 return self._basetype(string, shape, coefs) 

775 

776 @convert_operands(sameShape=True) 

777 @refine_operands(stop_at_affine=True) 

778 def __add__(self, other): 

779 """Denote addition from the right hand side.""" 

780 if not isinstance(other, BiaffineExpression): 

781 return NotImplemented 

782 

783 string = glyphs.clever_add(self.string, other.string) 

784 

785 coefs = {} 

786 for mtbs, coef in self._coefs.items(): 

787 coefs[mtbs] = coef + other._coefs[mtbs] \ 

788 if mtbs in other._coefs else coef 

789 for mtbs, coef in other._coefs.items(): 

790 coefs.setdefault(mtbs, coef) 

791 

792 return self._common_basetype(other)(string, self._shape, coefs) 

793 

794 @convert_operands(sameShape=True) 

795 @refine_operands(stop_at_affine=True) 

796 def __radd__(self, other): 

797 """Denote addition from the left hand side.""" 

798 if not isinstance(other, BiaffineExpression): 

799 return NotImplemented 

800 

801 return other.__add__(self) 

802 

803 @convert_operands(sameShape=True) 

804 @refine_operands(stop_at_affine=True) 

805 def __sub__(self, other): 

806 """Denote substraction from the right hand side.""" 

807 if not isinstance(other, BiaffineExpression): 

808 return NotImplemented 

809 

810 string = glyphs.clever_sub(self.string, other.string) 

811 

812 coefs = {} 

813 for mtbs, coef in self._coefs.items(): 

814 coefs[mtbs] = coef - other._coefs[mtbs] \ 

815 if mtbs in other._coefs else coef 

816 for mtbs, coef in other._coefs.items(): 

817 coefs.setdefault(mtbs, -coef) 

818 

819 return self._common_basetype(other)(string, self._shape, coefs) 

820 

821 @convert_operands(sameShape=True) 

822 @refine_operands(stop_at_affine=True) 

823 def __rsub__(self, other): 

824 """Denote substraction with self on the right hand side.""" 

825 if not isinstance(other, BiaffineExpression): 

826 return NotImplemented 

827 

828 return other.__sub__(self) 

829 

830 def __pos__(self): 

831 """Denote identity.""" 

832 return self 

833 

834 @cached_selfinverse_unary_operator 

835 def __neg__(self): 

836 """Denote negation.""" 

837 string = glyphs.clever_neg(self.string) 

838 coefs = {mtbs: -coef for mtbs, coef in self._coefs.items()} 

839 

840 return self._basetype(string, self._shape, coefs) 

841 

842 @convert_operands(rMatMul=True) 

843 @refine_operands(stop_at_affine=True) 

844 def __mul__(self, other): 

845 """Denote matrix multiplication from the right hand side.""" 

846 if not isinstance(other, BiaffineExpression): 

847 return NotImplemented 

848 

849 string = glyphs.clever_mul(self.string, other.string) 

850 

851 m, n = self._shape 

852 p, q = other._shape 

853 

854 # Handle or prepare multiplication with a scalar. 

855 if 1 in (m*n, p*q): 

856 if self.constant or other.constant: 

857 # Handle all cases involving a constant operand immediately. 

858 if self.constant: 

859 lhs = self._constant_coef 

860 RHS = other._coefs 

861 else: 

862 lhs = other._constant_coef 

863 RHS = self._coefs 

864 

865 shape = other._shape if m*n == 1 else self._shape 

866 

867 # Work around CVXOPT not broadcasting sparse scalars. 

868 if lhs.size == (1, 1): 

869 lhs = lhs[0] 

870 

871 return self._common_basetype(other)( 

872 string, shape, {mtbs: lhs*rhs for mtbs, rhs in RHS.items()}) 

873 elif n == p: 

874 # Regular matrix multiplication already works. 

875 pass 

876 elif m*n == 1: 

877 # Prepare a regular matrix multiplication. 

878 # HACK: Replaces self. 

879 self = self*cvxopt.spdiag([1]*p) 

880 m, n = self._shape 

881 else: 

882 # Prepare a regular matrix multiplication. 

883 assert p*q == 1 

884 other = other*cvxopt.spdiag([1]*n) 

885 p, q = other._shape 

886 

887 assert n == p 

888 

889 # Handle the common row by column multiplication more efficiently. 

890 # This includes some scalar by scalar products (both sides nonconstant). 

891 row_by_column = m == 1 and q == 1 

892 

893 # Regular matrix multiplication. 

894 coefs = {} 

895 for (lmtbs, lcoef), (rmtbs, rcoef) in product( 

896 self._coefs.items(), other._coefs.items()): 

897 if len(lmtbs) + len(rmtbs) > 2: 

898 raise NotImplementedError("Product not supported: " 

899 "One operand is biaffine and the other nonconstant.") 

900 elif len(lmtbs) == 0: 

901 # Constant by constant, linear or bilinear. 

902 if row_by_column: 

903 factor = lcoef.T 

904 else: 

905 # Use identity vec(AB) = (I ⊗ A)vec(B). 

906 factor = left_kronecker_I(lcoef, q, reshape=(m, n)) 

907 

908 self._update_coefs(coefs, rmtbs, factor*rcoef) 

909 elif len(rmtbs) == 0: 

910 # Constant, linear or bilinear by constant. 

911 if row_by_column: 

912 factor = rcoef.T 

913 else: 

914 # Use identity vec(AB) = (Bᵀ ⊗ I)vec(A). 

915 factor = right_kronecker_I( 

916 rcoef, m, reshape=(p, q), postT=True) 

917 

918 self._update_coefs(coefs, lmtbs, factor*lcoef) 

919 else: 

920 # Linear by linear. 

921 assert len(lmtbs) == 1 and len(rmtbs) == 1 

922 

923 if row_by_column: 

924 # Use identity (Ax)ᵀ(By) = vec(AᵀB)ᵀ·vec(xyᵀ). 

925 coef = (lcoef.T*rcoef)[:].T 

926 else: 

927 # Recipe found in "Robust conic optimization in Python" 

928 # (Stahlberg 2020). 

929 a, b = lmtbs[0].dim, rmtbs[0].dim 

930 

931 A = cvxopt_K(m, n)*lcoef 

932 B = rcoef 

933 

934 A = A.T 

935 B = B.T 

936 A.size = m*n*a, 1 

937 B.size = p*q*b, 1 

938 

939 A = cvxopt.spdiag([cvxopt_K(a, n)]*m)*A 

940 B = cvxopt.spdiag([cvxopt_K(b, p)]*q)*B 

941 A.size = n, m*a 

942 B.size = p, q*b 

943 A = A.T 

944 

945 C = A*B 

946 C.size = a, m*q*b 

947 C = C*cvxopt.spdiag([cvxopt_K(b, m)]*q) 

948 C.size = a*b, m*q 

949 C = C.T 

950 

951 coef = C 

952 

953 # Forbid quadratic results. 

954 if coef and lmtbs[0] is rmtbs[0]: 

955 raise TypeError("The product of {} and {} is quadratic " 

956 "in {} but a biaffine result is required here." 

957 .format(self.string, other.string, lmtbs[0].name)) 

958 

959 self._update_coefs(coefs, lmtbs + rmtbs, coef) 

960 

961 return self._common_basetype(other)(string, (m, q), coefs) 

962 

963 @convert_operands(lMatMul=True) 

964 @refine_operands(stop_at_affine=True) 

965 def __rmul__(self, other): 

966 """Denote matrix multiplication from the left hand side.""" 

967 if not isinstance(other, BiaffineExpression): 

968 return NotImplemented 

969 

970 return other.__mul__(self) 

971 

972 @convert_operands(sameShape=True) 

973 @refine_operands(stop_at_affine=True) 

974 def __or__(self, other): 

975 r"""Denote the scalar product with self on the left hand side. 

976 

977 For (complex) vectors :math:`a` and :math:`b` this is the dot product 

978 

979 .. math:: 

980 (a \mid b) 

981 &= \langle a, b \rangle \\ 

982 &= a \cdot b \\ 

983 &= b^H a. 

984 

985 For (complex) matrices :math:`A` and :math:`B` this is the Frobenius 

986 inner product 

987 

988 .. math:: 

989 (A \mid B) 

990 &= \langle A, B \rangle_F \\ 

991 &= A : B \\ 

992 &= \operatorname{tr}(B^H A) \\ 

993 &= \operatorname{vec}(B)^H \operatorname{vec}(\overline{A}) 

994 

995 .. note:: 

996 Write ``(A|B)`` instead of ``A|B`` for the scalar product of ``A`` 

997 and ``B`` to obtain correct operator binding within a larger 

998 expression context. 

999 """ 

1000 if not isinstance(other, BiaffineExpression): 

1001 return NotImplemented 

1002 

1003 result = other.vec.H * self.vec 

1004 result._symbStr = glyphs.clever_dotp( 

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

1006 return result 

1007 

1008 @convert_operands(sameShape=True) 

1009 @refine_operands(stop_at_affine=True) 

1010 def __ror__(self, other): 

1011 """Denote the scalar product with self on the right hand side.""" 

1012 if not isinstance(other, BiaffineExpression): 

1013 return NotImplemented 

1014 

1015 return other.__or__(self) 

1016 

1017 @convert_operands(sameShape=True) 

1018 @refine_operands(stop_at_affine=True) 

1019 def __xor__(self, other): 

1020 """Denote the elementwise (or Hadamard) product.""" 

1021 if not isinstance(other, BiaffineExpression): 

1022 return NotImplemented 

1023 

1024 string = glyphs.hadamard(self.string, other.string) 

1025 

1026 if other.constant: 

1027 factor = cvxopt.spdiag(other._constant_coef) 

1028 coefs = {mtbs: factor*coef for mtbs, coef in self._coefs.items()} 

1029 elif self.constant: 

1030 factor = cvxopt.spdiag(self._constant_coef) 

1031 coefs = {mtbs: factor*coef for mtbs, coef in other._coefs.items()} 

1032 else: 

1033 return (self.diag*other.vec).reshaped(self._shape).renamed(string) 

1034 

1035 return self._common_basetype(other)(string, self._shape, coefs) 

1036 

1037 @convert_operands(sameShape=True) 

1038 @refine_operands(stop_at_affine=True) 

1039 def __rxor__(self, other): 

1040 """Denote the elementwise (or Hadamard) product.""" 

1041 if not isinstance(other, BiaffineExpression): 

1042 return NotImplemented 

1043 

1044 return other.__xor__(self) 

1045 

1046 @convert_operands() 

1047 @refine_operands(stop_at_affine=True) 

1048 def __matmul__(self, other): 

1049 """Denote the Kronecker product from the right hand side.""" 

1050 if not isinstance(other, BiaffineExpression): 

1051 return NotImplemented 

1052 

1053 cls = self._common_basetype(other) 

1054 tc = cls._get_typecode() 

1055 

1056 string = glyphs.kron(self.string, other.string) 

1057 

1058 m, n = self._shape 

1059 p, q = other._shape 

1060 

1061 # Recipe in "Robust conic optimization in Python" (Stahlberg 2020). 

1062 Kqm = cvxopt_K(q, m) 

1063 Kqm_Ip = right_kronecker_I(Kqm, p) 

1064 In_Kqm_Ip = left_kronecker_I(Kqm_Ip, n) 

1065 

1066 def _kron(A, B): 

1067 A_B = load_data(numpy.kron(cvx2np(A), cvx2np(B)), typecode=tc)[0] 

1068 Kab = cvxopt_K(A.size[1], B.size[1]) 

1069 return In_Kqm_Ip*A_B*Kab 

1070 

1071 coefs = {} 

1072 for (lmtbs, lcoef), (rmtbs, rcoef) in product( 

1073 self._coefs.items(), other._coefs.items()): 

1074 if len(lmtbs) + len(rmtbs) <= 2: 

1075 if len(lmtbs) == 1 and len(rmtbs) == 1 and lmtbs[0] is rmtbs[0]: 

1076 raise TypeError("The product of {} and {} is quadratic " 

1077 "in {} but a biaffine result is required here." 

1078 .format(self.string, other.string, lmtbs[0].name)) 

1079 

1080 self._update_coefs(coefs, lmtbs + rmtbs, _kron(lcoef, rcoef)) 

1081 else: 

1082 raise NotImplementedError("Product not supported: " 

1083 "One operand is biaffine and the other nonconstant.") 

1084 

1085 return cls(string, (m*p, n*q), coefs) 

1086 

1087 @convert_operands() 

1088 @refine_operands(stop_at_affine=True) 

1089 def __rmatmul__(self, other): 

1090 """Denote the Kronecker product from the left hand side.""" 

1091 if not isinstance(other, BiaffineExpression): 

1092 return NotImplemented 

1093 

1094 return other.__matmul__(self) 

1095 

1096 @deprecated("2.2", "Use infix @ instead.") 

1097 def kron(self, other): 

1098 """Denote the Kronecker product from the right hand side.""" 

1099 return self.__matmul__(other) 

1100 

1101 @deprecated("2.2", "Reverse operands and use infix @ instead.") 

1102 def leftkron(self, other): 

1103 """Denote the Kronecker product from the left hand side.""" 

1104 return self.__rmatmul__(other) 

1105 

1106 @convert_operands(scalarRHS=True) 

1107 @refine_operands(stop_at_affine=True) 

1108 def __truediv__(self, other): 

1109 """Denote elementwise division by a scalar constant.""" 

1110 if not isinstance(other, BiaffineExpression): 

1111 return NotImplemented 

1112 

1113 if not other.constant: 

1114 raise TypeError( 

1115 "You may only divide {} by a constant.".format(self.string)) 

1116 

1117 if other.is0: 

1118 raise ZeroDivisionError( 

1119 "Tried to divide {} by zero.".format(self.string)) 

1120 

1121 divisor = other._constant_coef[0] 

1122 

1123 string = glyphs.div(self.string, other.string) 

1124 coefs = {mtbs: coef / divisor for mtbs, coef in self._coefs.items()} 

1125 

1126 return self._common_basetype(other)(string, self._shape, coefs) 

1127 

1128 @convert_operands(scalarLHS=True) 

1129 @refine_operands(stop_at_affine=True) 

1130 def __rtruediv__(self, other): 

1131 """Denote elementwise division by scalar constant self.""" 

1132 if not isinstance(other, BiaffineExpression): 

1133 return NotImplemented 

1134 

1135 return other.__truediv__(self) 

1136 

1137 @convert_operands(horiCat=True) 

1138 @refine_operands(stop_at_affine=True) 

1139 def __and__(self, other): 

1140 """Denote horizontal concatenation from the right hand side.""" 

1141 if not isinstance(other, BiaffineExpression): 

1142 return NotImplemented 

1143 

1144 string = glyphs.matrix_cat(self.string, other.string, horizontal=True) 

1145 shape = (self._shape[0], self._shape[1] + other._shape[1]) 

1146 

1147 coefs = {} 

1148 for mtbs in set(self._coefs.keys()).union(other._coefs.keys()): 

1149 l = self._coefs[mtbs] if mtbs in self._coefs else cvxopt.spmatrix( 

1150 [], [], [], (len(self), other._coefs[mtbs].size[1])) 

1151 r = other._coefs[mtbs] if mtbs in other._coefs else cvxopt.spmatrix( 

1152 [], [], [], (len(other), self._coefs[mtbs].size[1])) 

1153 

1154 coefs[mtbs] = cvxopt_vcat([l, r]) 

1155 

1156 return self._common_basetype(other)(string, shape, coefs) 

1157 

1158 @convert_operands(horiCat=True) 

1159 @refine_operands(stop_at_affine=True) 

1160 def __rand__(self, other): 

1161 """Denote horizontal concatenation from the left hand side.""" 

1162 if not isinstance(other, BiaffineExpression): 

1163 return NotImplemented 

1164 

1165 return other.__and__(self) 

1166 

1167 @convert_operands(vertCat=True) 

1168 @refine_operands(stop_at_affine=True) 

1169 def __floordiv__(self, other): 

1170 """Denote vertical concatenation from below.""" 

1171 def interleave_columns(upper, lower, upperRows, lowerRows, cols): 

1172 p, q = upperRows, lowerRows 

1173 return [column for columnPairs in [ 

1174 (upper[j*p:j*p+p, :], lower[j*q:j*q+q, :]) for j in range(cols)] 

1175 for column in columnPairs] 

1176 

1177 if not isinstance(other, BiaffineExpression): 

1178 return NotImplemented 

1179 

1180 string = glyphs.matrix_cat(self.string, other.string, horizontal=False) 

1181 

1182 p, q, c = self._shape[0], other._shape[0], self._shape[1] 

1183 shape = (p + q, c) 

1184 

1185 coefs = {} 

1186 for mtbs in set(self._coefs.keys()).union(other._coefs.keys()): 

1187 u = self._coefs[mtbs] if mtbs in self._coefs else cvxopt.spmatrix( 

1188 [], [], [], (len(self), other._coefs[mtbs].size[1])) 

1189 l = other._coefs[mtbs] if mtbs in other._coefs else cvxopt.spmatrix( 

1190 [], [], [], (len(other), self._coefs[mtbs].size[1])) 

1191 

1192 coefs[mtbs] = cvxopt_vcat(interleave_columns(u, l, p, q, c)) 

1193 

1194 return self._common_basetype(other)(string, shape, coefs) 

1195 

1196 @convert_operands(vertCat=True) 

1197 @refine_operands(stop_at_affine=True) 

1198 def __rfloordiv__(self, other): 

1199 """Denote vertical concatenation from above.""" 

1200 if not isinstance(other, BiaffineExpression): 

1201 return NotImplemented 

1202 

1203 return other.__floordiv__(self) 

1204 

1205 @convert_operands(scalarRHS=True) 

1206 @refine_operands() # Refine both sides to real if possible. 

1207 def __pow__(self, other): 

1208 """Denote exponentiation.""" 

1209 from .exp_powtrace import PowerTrace 

1210 

1211 if not self.scalar: 

1212 raise TypeError("May only exponentiate a scalar expression.") 

1213 

1214 if not other.constant: 

1215 raise TypeError("The exponent must be constant.") 

1216 

1217 if other.complex: 

1218 raise TypeError("The exponent must be real-valued.") 

1219 

1220 exponent = other.value 

1221 

1222 if exponent == 2: 

1223 return (self | self) # Works for complex base. 

1224 else: 

1225 return PowerTrace(self, exponent) # Errors on complex base. 

1226 

1227 # -------------------------------------------------------------------------- 

1228 # Properties and functions that describe the expression. 

1229 # -------------------------------------------------------------------------- 

1230 

1231 @cached_property 

1232 def hermitian(self): # noqa (D402 thinks this includes a signature) 

1233 """Whether the expression is a hermitian (or symmetric) matrix. 

1234 

1235 Uses :data:`~picos.settings.RELATIVE_HERMITIANNESS_TOLERANCE`. 

1236 

1237 If PICOS rejects your near-hermitian (near-symmetric) expression as not 

1238 hermitian (not symmetric), you can use :meth:`hermitianized` to correct 

1239 larger numeric errors or the effects of noisy data. 

1240 """ 

1241 return self.equals( 

1242 self.H, relTol=settings.RELATIVE_HERMITIANNESS_TOLERANCE) 

1243 

1244 @property 

1245 def is0(self): 

1246 """Whether this is a constant scalar, vector or matrix of all zeros.""" 

1247 return not self._coefs 

1248 

1249 @cached_property 

1250 def is1(self): 

1251 """Whether this is a constant scalar or vector of all ones.""" 

1252 # Must be a scalar or vector. 

1253 if self._shape[0] != 1 and self._shape[1] != 1: 

1254 return False 

1255 

1256 # Must be constant. 

1257 if not self.constant: 

1258 return False 

1259 

1260 # Constant term must be all ones. 

1261 return not self._constant_coef - 1 

1262 

1263 @cached_property 

1264 def isI(self): 

1265 """Whether this is a constant identity matrix.""" 

1266 m, n = self._shape 

1267 

1268 # Must be a square matrix. 

1269 if m != n: 

1270 return False 

1271 

1272 # Must be constant. 

1273 if not self.constant: 

1274 return False 

1275 

1276 # Constant term must be the identity. 

1277 return not self._constant_coef - cvxopt.spdiag([1]*m)[:] 

1278 

1279 @cached_property 

1280 def complex(self): 

1281 """Whether the expression can be complex-valued.""" 

1282 # NOTE: The typecode check works around a bug in CVXOPT prior to 1.2.4. 

1283 return any(coef.typecode == "z" and coef.imag() 

1284 for coef in self._coefs.values()) 

1285 

1286 @property 

1287 def isreal(self): 

1288 """Whether the expression is always real-valued.""" 

1289 return not self.complex 

1290 

1291 @convert_operands() 

1292 def equals(self, other, absTol=None, relTol=None): 

1293 """Check mathematical equality with another (bi)affine expression. 

1294 

1295 The precise type of both (bi)affine expressions may differ. In 

1296 particular, a :class:`~.exp_affine.ComplexAffineExpression` with real 

1297 coefficients and constant term can be equal to an 

1298 :class:`~.exp_affine.AffineExpression`. 

1299 

1300 If the operand is not already a PICOS expression, an attempt is made to 

1301 load it as a constant affine expression. In this case, no reshaping or 

1302 broadcasting is used to bring the constant term to the same shape as 

1303 this expression. In particular, 

1304 

1305 - ``0`` refers to a scalar zero (see also :meth:`is0`), 

1306 - lists and tuples are treated as column vectors and 

1307 - algebraic strings must specify a shape (see 

1308 :func:`~.data.load_data`). 

1309 

1310 :param other: Another PICOS expression or a constant numeric data value 

1311 supported by :func:`~.data.load_data`. 

1312 :param absTol: As long as all absolute differences between scalar 

1313 entries of the coefficient matrices and the constant terms being 

1314 compared does not exceed this bound, consider the expressions equal. 

1315 :param relTol: As long as all absolute differences between scalar 

1316 entries of the coefficient matrices and the constant terms being 

1317 compared divided by the maximum absolute value found in either term 

1318 does not exceed this bound, consider the expressions equal. 

1319 

1320 :Example: 

1321 

1322 >>> from picos import Constant 

1323 >>> A = Constant("A", 0, (5,5)) 

1324 >>> repr(A) 

1325 '<5×5 Real Constant: A>' 

1326 >>> A.is0 

1327 True 

1328 >>> A.equals(0) 

1329 False 

1330 >>> A.equals("|0|(5,5)") 

1331 True 

1332 >>> repr(A*1j) 

1333 '<5×5 Complex Constant: A·1j>' 

1334 >>> A.equals(A*1j) 

1335 True 

1336 """ 

1337 if self is other: 

1338 return True 

1339 

1340 if not isinstance(other, BiaffineExpression): 

1341 return False 

1342 

1343 if self._shape != other._shape: 

1344 return False 

1345 

1346 assert not any((mtbs[1], mtbs[0]) in other._bilinear_coefs 

1347 for mtbs in self._bilinear_coefs), \ 

1348 "{} and {} store bilinear terms in a different order." \ 

1349 .format(self.string, other.string) 

1350 

1351 for mtbs in other._coefs: 

1352 if mtbs not in self._coefs: 

1353 return False 

1354 

1355 for mtbs, coef in self._coefs.items(): 

1356 if mtbs not in other._coefs: 

1357 return False 

1358 

1359 if not cvxopt_equals(coef, other._coefs[mtbs], absTol, relTol): 

1360 return False 

1361 

1362 return True 

1363 

1364 # -------------------------------------------------------------------------- 

1365 # Methods and properties that return modified copies. 

1366 # -------------------------------------------------------------------------- 

1367 

1368 def renamed(self, string): 

1369 """Return the expression with a modified string description.""" 

1370 return self._basetype(string, self._shape, self._coefs) 

1371 

1372 def reshaped(self, shape, order="F"): 

1373 r"""Return the expression reshaped in the given order. 

1374 

1375 The default indexing order is column-major. Given an :math:`m \times n` 

1376 matrix, reshaping in the default order is a constant time operation 

1377 while reshaping in row-major order requires :math:`O(mn)` time. However, 

1378 the latter allows you to be consistent with NumPy, which uses C-order (a 

1379 generalization of row-major) by default. 

1380 

1381 :param str order: 

1382 The indexing order to use when reshaping. Must be either ``"F"`` for 

1383 Fortran-order (column-major) or ``"C"`` for C-order (row-major). 

1384 

1385 :Example: 

1386 

1387 >>> from picos import Constant 

1388 >>> C = Constant("C", range(6), (2, 3)) 

1389 >>> print(C) 

1390 [ 0.00e+00 2.00e+00 4.00e+00] 

1391 [ 1.00e+00 3.00e+00 5.00e+00] 

1392 >>> print(C.reshaped((3, 2))) 

1393 [ 0.00e+00 3.00e+00] 

1394 [ 1.00e+00 4.00e+00] 

1395 [ 2.00e+00 5.00e+00] 

1396 >>> print(C.reshaped((3, 2), order="C")) 

1397 [ 0.00e+00 2.00e+00] 

1398 [ 4.00e+00 1.00e+00] 

1399 [ 3.00e+00 5.00e+00] 

1400 """ 

1401 if order not in "FC": 

1402 raise ValueError("Order must be given as 'F' or 'C'.") 

1403 

1404 shape = load_shape(shape, wildcards=True) 

1405 

1406 if shape == self._shape: 

1407 return self 

1408 elif shape == (None, None): 

1409 return self 

1410 

1411 length = len(self) 

1412 

1413 if shape[0] is None: 

1414 shape = (length // shape[1], shape[1]) 

1415 elif shape[1] is None: 

1416 shape = (shape[0], length // shape[0]) 

1417 

1418 if shape[0]*shape[1] != length: 

1419 raise ValueError("Can only reshape to a matrix of same size.") 

1420 

1421 if order == "F": 

1422 coefs = self._coefs 

1423 reshaped_glyph = glyphs.reshaped 

1424 else: 

1425 m, n = self._shape 

1426 p, q = shape 

1427 

1428 K_old = cvxopt_K(m, n, self._typecode) 

1429 K_new = cvxopt_K(q, p, self._typecode) 

1430 R = K_new * K_old 

1431 

1432 coefs = {mtbs: R * coef for mtbs, coef in self._coefs.items()} 

1433 reshaped_glyph = glyphs.reshaprm 

1434 

1435 string = reshaped_glyph(self.string, glyphs.shape(shape)) 

1436 return self._basetype(string, shape, coefs) 

1437 

1438 def broadcasted(self, shape): 

1439 """Return the expression broadcasted to the given shape. 

1440 

1441 :Example: 

1442 

1443 >>> from picos import Constant 

1444 >>> C = Constant("C", range(6), (2, 3)) 

1445 >>> print(C) 

1446 [ 0.00e+00 2.00e+00 4.00e+00] 

1447 [ 1.00e+00 3.00e+00 5.00e+00] 

1448 >>> print(C.broadcasted((6, 6))) 

1449 [ 0.00e+00 2.00e+00 4.00e+00 0.00e+00 2.00e+00 4.00e+00] 

1450 [ 1.00e+00 3.00e+00 5.00e+00 1.00e+00 3.00e+00 5.00e+00] 

1451 [ 0.00e+00 2.00e+00 4.00e+00 0.00e+00 2.00e+00 4.00e+00] 

1452 [ 1.00e+00 3.00e+00 5.00e+00 1.00e+00 3.00e+00 5.00e+00] 

1453 [ 0.00e+00 2.00e+00 4.00e+00 0.00e+00 2.00e+00 4.00e+00] 

1454 [ 1.00e+00 3.00e+00 5.00e+00 1.00e+00 3.00e+00 5.00e+00] 

1455 """ 

1456 shape = load_shape(shape, wildcards=True) 

1457 shape = blend_shapes(shape, self._shape) 

1458 

1459 if shape == self._shape: 

1460 return self 

1461 

1462 vdup = shape[0] // self._shape[0] 

1463 hdup = shape[1] // self._shape[1] 

1464 

1465 if (self._shape[0] * vdup, self._shape[1] * hdup) != shape: 

1466 raise ValueError("Cannot broadcast from shape {} to {}." 

1467 .format(glyphs.shape(self._shape), glyphs.shape(shape))) 

1468 

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

1470 string = glyphs.matrix(self.string) 

1471 return (self * cvxopt.matrix(1.0, shape)).renamed(string) 

1472 

1473 string = glyphs.bcasted(self.string, glyphs.shape(shape)) 

1474 return (cvxopt.matrix(1.0, (vdup, hdup)) @ self).renamed(string) 

1475 

1476 def reshaped_or_broadcasted(self, shape): 

1477 """Return the expression :meth:`reshaped` or :meth:`broadcasted`. 

1478 

1479 Unlike with :meth:`reshaped` and :meth:`broadcasted`, the target shape 

1480 may not contain a wildcard character. 

1481 

1482 If the wildcard-free target shape has the same number of elements as 

1483 the current shape, then this is the same as :meth:`reshaped`, otherwise 

1484 it is the same as :meth:`broadcasted`. 

1485 """ 

1486 shape = load_shape(shape, wildcards=False) 

1487 

1488 try: 

1489 if shape[0]*shape[1] == len(self): 

1490 return self.reshaped(shape) 

1491 else: 

1492 return self.broadcasted(shape) 

1493 except ValueError: 

1494 raise ValueError( 

1495 "Cannot reshape or broadcast from shape {} to {}.".format( 

1496 glyphs.shape(self._shape), glyphs.shape(shape))) from None 

1497 

1498 @cached_property 

1499 def hermitianized(self): 

1500 r"""The expression projected onto the subspace of hermitian matrices. 

1501 

1502 For a square (complex) affine expression :math:`A`, this is 

1503 :math:`\frac{1}{2}(A + A^H)`. 

1504 

1505 If the expression is not complex, then this is a projection onto the 

1506 subspace of symmetric matrices. 

1507 """ 

1508 if not self.square: 

1509 raise TypeError("Cannot hermitianize non-square {}.".format(self)) 

1510 

1511 return (self + self.H)/2 

1512 

1513 @cached_property 

1514 def real(self): 

1515 """Real part of the expression.""" 

1516 return self._basetype(glyphs.real(self.string), self._shape, 

1517 {mtbs: coef.real() for mtbs, coef in self._coefs.items()}) 

1518 

1519 @cached_property 

1520 def imag(self): 

1521 """Imaginary part of the expression.""" 

1522 return self._basetype(glyphs.imag(self.string), self._shape, 

1523 {mtbs: coef.imag() for mtbs, coef in self._coefs.items()}) 

1524 

1525 @cached_property 

1526 def bilin(self): 

1527 """Pure bilinear part of the expression.""" 

1528 return self._basetype( 

1529 glyphs.blinpart(self._symbStr), self._shape, self._bilinear_coefs) 

1530 

1531 @cached_property 

1532 def lin(self): 

1533 """Linear part of the expression.""" 

1534 return self._basetype( 

1535 glyphs.linpart(self._symbStr), self._shape, self._linear_coefs) 

1536 

1537 @cached_property 

1538 def noncst(self): 

1539 """Nonconstant part of the expression.""" 

1540 coefs = {mtbs: coefs for mtbs, coefs in self._coefs.items() if mtbs} 

1541 

1542 return self._basetype( 

1543 glyphs.ncstpart(self._symbStr), self._shape, coefs) 

1544 

1545 @cached_property 

1546 def cst(self): 

1547 """Constant part of the expression.""" 

1548 coefs = {(): self._coefs[()]} if () in self._coefs else {} 

1549 

1550 return self._basetype(glyphs.cstpart(self._symbStr), self._shape, coefs) 

1551 

1552 @cached_selfinverse_property 

1553 def T(self): 

1554 """Matrix transpose.""" 

1555 if len(self) == 1: 

1556 return self 

1557 

1558 m, n = self._shape 

1559 K = cvxopt_K(m, n, self._typecode) 

1560 

1561 string = glyphs.transp(self.string) 

1562 shape = (self._shape[1], self._shape[0]) 

1563 coefs = {mtbs: K * coef for mtbs, coef in self._coefs.items()} 

1564 

1565 return self._basetype(string, shape, coefs) 

1566 

1567 @cached_selfinverse_property 

1568 def conj(self): 

1569 """Complex conjugate.""" 

1570 string = glyphs.conj(self.string) 

1571 coefs = {mtbs: coef.H.T for mtbs, coef in self._coefs.items()} 

1572 

1573 return self._basetype(string, self._shape, coefs) 

1574 

1575 @cached_selfinverse_property 

1576 def H(self): 

1577 """Conjugate (or Hermitian) transpose.""" 

1578 return self.T.conj.renamed(glyphs.htransp(self._symbStr)) 

1579 

1580 def _square_equal_subsystem_dims(self, diagLen): 

1581 """Support :func:`partial_trace` and :func:`partial_transpose`.""" 

1582 m, n = self._shape 

1583 k = math.log(m, diagLen) 

1584 

1585 if m != n or int(k) != k: 

1586 raise TypeError("The expression has shape {} so it cannot be " 

1587 "decomposed into subsystems of shape {}.".format( 

1588 glyphs.shape(self._shape), glyphs.shape((diagLen,)*2))) 

1589 

1590 return ((diagLen,)*2,)*int(k) 

1591 

1592 def partial_transpose(self, subsystems, dimensions=2): 

1593 r"""Return the expression with selected subsystems transposed. 

1594 

1595 If the expression can be written as 

1596 :math:`A_0 \otimes \cdots \otimes A_{n-1}` for matrices 

1597 :math:`A_0, \ldots, A_{n-1}` with shapes given in ``dimensions``, then 

1598 this returns :math:`B_0 \otimes \cdots \otimes B_{n-1}` with 

1599 :math:`B_i = A_i^T`, if ``i in subsystems`` (with :math:`i = -1` read as 

1600 :math:`n-1`), and :math:`B_i = A_i`, otherwise. 

1601 

1602 :param subsystems: A collection of or a single subystem number, indexed 

1603 from zero, corresponding to subsystems that shall be transposed. 

1604 The value :math:`-1` refers to the last subsystem. 

1605 :type subsystems: int or tuple or list 

1606 

1607 :param dimensions: Either an integer :math:`d` so that the subsystems 

1608 are assumed to be all of shape :math:`d \times d`, or a sequence of 

1609 subsystem shapes where an integer :math:`d` within the sequence is 

1610 read as :math:`d \times d`. In any case, the elementwise product 

1611 over all subsystem shapes must equal the expression's shape. 

1612 :type dimensions: int or tuple or list 

1613 

1614 :raises TypeError: If the subsystems do not match the expression. 

1615 :raises IndexError: If the subsystem selection is invalid. 

1616 

1617 :Example: 

1618 

1619 >>> from picos import Constant 

1620 >>> A = Constant("A", range(16), (4, 4)) 

1621 >>> print(A) #doctest: +NORMALIZE_WHITESPACE 

1622 [ 0.00e+00 4.00e+00 8.00e+00 1.20e+01] 

1623 [ 1.00e+00 5.00e+00 9.00e+00 1.30e+01] 

1624 [ 2.00e+00 6.00e+00 1.00e+01 1.40e+01] 

1625 [ 3.00e+00 7.00e+00 1.10e+01 1.50e+01] 

1626 >>> A0 = A.partial_transpose(0); A0 

1627 <4×4 Real Constant: A.{[2×2]ᵀ⊗[2×2]}> 

1628 >>> print(A0) #doctest: +NORMALIZE_WHITESPACE 

1629 [ 0.00e+00 4.00e+00 2.00e+00 6.00e+00] 

1630 [ 1.00e+00 5.00e+00 3.00e+00 7.00e+00] 

1631 [ 8.00e+00 1.20e+01 1.00e+01 1.40e+01] 

1632 [ 9.00e+00 1.30e+01 1.10e+01 1.50e+01] 

1633 >>> A1 = A.partial_transpose(1); A1 

1634 <4×4 Real Constant: A.{[2×2]⊗[2×2]ᵀ}> 

1635 >>> print(A1) #doctest: +NORMALIZE_WHITESPACE 

1636 [ 0.00e+00 1.00e+00 8.00e+00 9.00e+00] 

1637 [ 4.00e+00 5.00e+00 1.20e+01 1.30e+01] 

1638 [ 2.00e+00 3.00e+00 1.00e+01 1.10e+01] 

1639 [ 6.00e+00 7.00e+00 1.40e+01 1.50e+01] 

1640 """ 

1641 m, n = self._shape 

1642 

1643 if isinstance(dimensions, int): 

1644 dimensions = self._square_equal_subsystem_dims(dimensions) 

1645 else: 

1646 dimensions = [ 

1647 (d, d) if isinstance(d, int) else d for d in dimensions] 

1648 

1649 if reduce( 

1650 lambda x, y: (x[0]*y[0], x[1]*y[1]), dimensions) != self._shape: 

1651 raise TypeError("Subsystem dimensions do not match expression.") 

1652 

1653 if isinstance(subsystems, int): 

1654 subsystems = (subsystems,) 

1655 elif not subsystems: 

1656 return self 

1657 

1658 numSys = len(dimensions) 

1659 subsystems = set(numSys - 1 if sys == -1 else sys for sys in subsystems) 

1660 

1661 for sys in subsystems: 

1662 if not isinstance(sys, int): 

1663 raise IndexError("Subsystem indices must be integer, not {}." 

1664 .format(type(sys).__name__)) 

1665 elif sys < 0: 

1666 raise IndexError("Subsystem indices must be nonnegative.") 

1667 elif sys >= numSys: 

1668 raise IndexError( 

1669 "Subsystem index {} out of range for {} systems total." 

1670 .format(sys, numSys)) 

1671 

1672 # If all subsystems are transposed, this is regular transposition. 

1673 if len(subsystems) == numSys: 

1674 return self.T 

1675 

1676 # Prepare sparse K such that K·vec(A) = vec(partial_transpose(A)). 

1677 d = m * n 

1678 V = [1]*d 

1679 I = range(d) 

1680 J = cvxopt.matrix(I) 

1681 T = cvxopt.matrix(0, J.size) 

1682 obh, obw = 1, 1 

1683 sysStrings = None 

1684 

1685 # Apply transpositions starting with the rightmost Kronecker factor. 

1686 for sys in range(numSys - 1, -1, -1): 

1687 # Shape of current system. 

1688 p, q = dimensions[sys] 

1689 sysString = glyphs.matrix(glyphs.shape((p, q))) 

1690 

1691 # Height/width of "inner" blocks being moved, initially scalars. 

1692 ibh, ibw = obh, obw 

1693 

1694 # Heigh/width of "outer" blocks whose relative position is 

1695 # maintained but that are subject to transposition independently. 

1696 # In the last iteration this is the shape of the resulting matrix. 

1697 obh *= p 

1698 obw *= q 

1699 

1700 # Only transpose selected subsystems. 

1701 if sys not in subsystems: 

1702 sysStrings = glyphs.kron(sysString, sysStrings) \ 

1703 if sysStrings else sysString 

1704 continue 

1705 else: 

1706 sysStrings = glyphs.kron(glyphs.transp(sysString), sysStrings) \ 

1707 if sysStrings else glyphs.transp(sysString) 

1708 

1709 # Shape of outer blocks after transposition. 

1710 obhT, obwT = obw // ibw * ibh, obh // ibh * ibw 

1711 

1712 # Shape of full matrix after transposition. 

1713 mT, nT = m // obh * obhT, n // obw * obwT 

1714 

1715 for vi in I: 

1716 # Full matrix column and row. 

1717 c, r = divmod(vi, m) 

1718 

1719 # Outer block vertical index and row within outer block, 

1720 # outer block horizontal index and column within outer block. 

1721 obi, obr = divmod(r, obh) 

1722 obj, obc = divmod(c, obw) 

1723 

1724 # Inner block vertical index and row within inner block, 

1725 # inner block horizontal index and column within inner block. 

1726 ibi, ibr = divmod(obr, ibh) 

1727 ibj, ibc = divmod(obc, ibw) 

1728 

1729 # (1) ibi*ibw + ibc is column within the transposed outer block; 

1730 # adding obj*obwT yields the column in the transposed matrix. 

1731 # (2) ibj*ibh + ibr is row within the transposed outer block; 

1732 # adding obi*obhT yields the row in the transposed matrix. 

1733 # (3) tvi is index within the vectorized transposed matrix. 

1734 tvi = (obj*obwT + ibi*ibw + ibc)*mT \ 

1735 + (obi*obhT + ibj*ibh + ibr) 

1736 

1737 # Prepare the transposition. 

1738 T[tvi] = J[vi] 

1739 

1740 # Apply the transposition. 

1741 J, T = T, J 

1742 m, n, obh, obw = mT, nT, obhT, obwT 

1743 

1744 # Finalize the partial commutation matrix. 

1745 K = cvxopt.spmatrix(V, I, J, (d, d), self._typecode) 

1746 

1747 string = glyphs.ptransp_(self.string, sysStrings) 

1748 shape = (m, n) 

1749 coefs = {mtbs: K * coef for mtbs, coef in self._coefs.items()} 

1750 

1751 return self._basetype(string, shape, coefs) 

1752 

1753 @cached_property 

1754 def T0(self): 

1755 r"""Expression with the first :math:`2 \times 2` subsystem transposed. 

1756 

1757 Only available for a :math:`2^k \times 2^k` matrix with all subsystems 

1758 of shape :math:`2 \times 2`. Use :meth:`partial_transpose` otherwise. 

1759 """ 

1760 return self.partial_transpose(subsystems=0) 

1761 

1762 @cached_property 

1763 def T1(self): 

1764 r"""Expression with the second :math:`2 \times 2` subsystem transposed. 

1765 

1766 Only available for a :math:`2^k \times 2^k` matrix with all subsystems 

1767 of shape :math:`2 \times 2`. Use :meth:`partial_transpose` otherwise. 

1768 """ 

1769 return self.partial_transpose(subsystems=1) 

1770 

1771 @cached_property 

1772 def T2(self): 

1773 r"""Expression with the third :math:`2 \times 2` subsystem transposed. 

1774 

1775 Only available for a :math:`2^k \times 2^k` matrix with all subsystems 

1776 of shape :math:`2 \times 2`. Use :meth:`partial_transpose` otherwise. 

1777 """ 

1778 return self.partial_transpose(subsystems=2) 

1779 

1780 @cached_property 

1781 def T3(self): 

1782 r"""Expression with the fourth :math:`2 \times 2` subsystem transposed. 

1783 

1784 Only available for a :math:`2^k \times 2^k` matrix with all subsystems 

1785 of shape :math:`2 \times 2`. Use :meth:`partial_transpose` otherwise. 

1786 """ 

1787 return self.partial_transpose(subsystems=3) 

1788 

1789 @cached_property 

1790 def Tl(self): 

1791 r"""Expression with the last :math:`2 \times 2` subsystem transposed. 

1792 

1793 Only available for a :math:`2^k \times 2^k` matrix with all subsystems 

1794 of shape :math:`2 \times 2`. Use :meth:`partial_transpose` otherwise. 

1795 """ 

1796 return self.partial_transpose(subsystems=-1) 

1797 

1798 @staticmethod 

1799 def _reindex_F(indices, source, destination): 

1800 """Convert indices between different tensor shapes in Fortran-order.""" 

1801 new = [] 

1802 offset = 0 

1803 factor = 1 

1804 

1805 for index, base in zip(indices, source): 

1806 offset += factor*index 

1807 factor *= base 

1808 

1809 for base in destination: 

1810 offset, remainder = divmod(offset, base) 

1811 new.append(remainder) 

1812 

1813 return tuple(new) 

1814 

1815 @staticmethod 

1816 def _reindex_C(indices, source, destination): 

1817 """Convert indices between different tensor shapes in C-order.""" 

1818 new = [] 

1819 offset = 0 

1820 factor = 1 

1821 

1822 for index, base in zip(reversed(indices), reversed(source)): 

1823 offset += factor*index 

1824 factor *= base 

1825 

1826 for base in reversed(destination): 

1827 offset, remainder = divmod(offset, base) 

1828 new.insert(0, remainder) 

1829 

1830 return tuple(new) 

1831 

1832 def reshuffled(self, permutation="ikjl", dimensions=None, order="C"): 

1833 """Return the reshuffled or realigned expression. 

1834 

1835 This operation works directly on matrices. However, it is equivalent to 

1836 the following sequence of operations: 

1837 

1838 1. The matrix is reshaped to a tensor with the given ``dimensions`` and 

1839 according to ``order``. 

1840 2. The tensor's axes are permuted according to ``permutation``. 

1841 3. The tensor is reshaped back to the shape of the original matrix 

1842 according to ``order``. 

1843 

1844 For comparison, the following function applies the same operation to a 

1845 2D NumPy :class:`~numpy:numpy.ndarray`: 

1846 

1847 .. code:: 

1848 

1849 def reshuffle_numpy(matrix, permutation, dimensions, order): 

1850 P = "{} -> {}".format("".join(sorted(permutation)), permutation) 

1851 reshuffled = numpy.reshape(matrix, dimensions, order) 

1852 reshuffled = numpy.einsum(P, reshuffled) 

1853 return numpy.reshape(reshuffled, matrix.shape, order) 

1854 

1855 :param permutation: 

1856 A sequence of comparable elements with length equal to the number of 

1857 tensor dimensions. The sequence is compared to its ordered version 

1858 and the resulting permutation pattern is used to permute the tensor 

1859 indices. For instance, the string ``"ikjl"`` is compared to its 

1860 sorted version ``"ijkl"`` and denotes that the second and third axis 

1861 should be swapped. 

1862 :type permutation: 

1863 str or tuple or list 

1864 

1865 :param dimensions: 

1866 If this is an integer sequence, then it defines the dimensions of 

1867 the tensor. If this is :obj:`None`, then the tensor is assumed to be 

1868 hypercubic and the number of dimensions is inferred from the 

1869 ``permutation`` argument. 

1870 :type dimensions: 

1871 None or tuple or list 

1872 

1873 :param str order: 

1874 The indexing order to use for the virtual reshaping. Must be either 

1875 ``"F"`` for Fortran-order (generalization of column-major) or 

1876 ``"C"`` for C-order (generalization of row-major). Note that PICOS 

1877 usually reshapes in Fortran-order while NumPy defaults to C-order. 

1878 

1879 :Example: 

1880 

1881 >>> from picos import Constant 

1882 >>> A = Constant("A", range(16), (4, 4)) 

1883 >>> print(A) #doctest: +NORMALIZE_WHITESPACE 

1884 [ 0.00e+00 4.00e+00 8.00e+00 1.20e+01] 

1885 [ 1.00e+00 5.00e+00 9.00e+00 1.30e+01] 

1886 [ 2.00e+00 6.00e+00 1.00e+01 1.40e+01] 

1887 [ 3.00e+00 7.00e+00 1.10e+01 1.50e+01] 

1888 >>> R = A.reshuffled(); R 

1889 <4×4 Real Constant: shuffled(A,ikjl,C)> 

1890 >>> print(R) #doctest: +NORMALIZE_WHITESPACE 

1891 [ 0.00e+00 4.00e+00 1.00e+00 5.00e+00] 

1892 [ 8.00e+00 1.20e+01 9.00e+00 1.30e+01] 

1893 [ 2.00e+00 6.00e+00 3.00e+00 7.00e+00] 

1894 [ 1.00e+01 1.40e+01 1.10e+01 1.50e+01] 

1895 >>> A.reshuffled("ji").equals(A.T) # Regular transposition. 

1896 True 

1897 >>> A.reshuffled("3214").equals(A.T0) # Partial transposition (1). 

1898 True 

1899 >>> A.reshuffled("1432").equals(A.T1) # Partial transposition (2). 

1900 True 

1901 """ 

1902 m, n = self._shape 

1903 mn = m*n 

1904 

1905 # Load the permutation. 

1906 ordered = sorted(permutation) 

1907 P = dict(enumerate(ordered.index(element) for element in permutation)) 

1908 

1909 if len(set(P.values())) < len(P): 

1910 raise ValueError("The sequence defining the permutation appears to " 

1911 "contain duplicate elements.") 

1912 

1913 assert not set(P.keys()).symmetric_difference(set(P.values())) 

1914 

1915 numDims = len(P) 

1916 

1917 # Load the dimensions. 

1918 guessDimensions = dimensions is None 

1919 

1920 if guessDimensions: 

1921 dimensions = (int(mn**(1.0 / numDims)),)*numDims 

1922 else: 

1923 if len(dimensions) != numDims: 

1924 raise ValueError("The number of indices does not match the " 

1925 "number of dimensions.") 

1926 

1927 if reduce(int.__mul__, dimensions, 1) != mn: 

1928 raise TypeError("The {} matrix {} cannot be reshaped to a {} " 

1929 "tensor.".format(glyphs.shape(self.shape), self.string, 

1930 "hypercubic order {}".format(numDims) if guessDimensions 

1931 else glyphs.size("", "").join(str(d) for d in dimensions))) 

1932 

1933 # Load the indexing order. 

1934 if order not in "FC": 

1935 raise ValueError("Order must be given as 'F' or 'C'.") 

1936 

1937 reindex = self._reindex_F if order == "F" else self._reindex_C 

1938 

1939 # Nothing to do for the neutral permutation. 

1940 if all(key == val for key, val in P.items()): 

1941 return self 

1942 

1943 # Create a sparse mtrix R such that R·vec(A) = vec(reshuffle(A)). 

1944 V, I, J = [1]*mn, [], range(mn) 

1945 for i in range(mn): 

1946 (k, j) = divmod(i, m) # (j, k) are column-major matrix indices. 

1947 indices = reindex((j, k), (m, n), dimensions) 

1948 newIndices = tuple(indices[P[d]] for d in range(numDims)) 

1949 newDimensions = tuple(dimensions[P[d]] for d in range(numDims)) 

1950 (j, k) = reindex(newIndices, newDimensions, (m, n)) 

1951 I.append(k*m + j) 

1952 R = cvxopt.spmatrix(V, I, J, (mn, mn), self._typecode) 

1953 

1954 # Create the string. 

1955 strArgs = [self.string, str(permutation).replace(" ", ""), order] 

1956 

1957 if not guessDimensions: 

1958 strArgs.insert(2, str(dimensions).replace(" ", "")) 

1959 

1960 string = glyphs.shuffled(",".join(strArgs)) 

1961 

1962 # Finalize the new expression. 

1963 shape = (m, n) 

1964 coefs = {mtbs: R * coef for mtbs, coef in self._coefs.items()} 

1965 

1966 return self._basetype(string, shape, coefs) 

1967 

1968 @cached_property 

1969 def sum(self): 

1970 """Sum over all scalar elements of the expression.""" 

1971 # NOTE: glyphs.clever_dotp detects this case and uses the sum glyph. 

1972 # NOTE: 1 on the right hand side in case self is complex. 

1973 return (self | 1) 

1974 

1975 @cached_property 

1976 def tr(self): 

1977 """Trace of a square expression.""" 

1978 if not self.square: 

1979 raise TypeError("Cannot compute the trace of non-square {}." 

1980 .format(self.string)) 

1981 

1982 # NOTE: glyphs.clever_dotp detects this case and uses the trace glyph. 

1983 # NOTE: "I" on the right hand side in case self is complex. 

1984 return (self | "I") 

1985 

1986 def partial_trace(self, subsystems, dimensions=2): 

1987 r"""Return the partial trace over selected subsystems. 

1988 

1989 If the expression can be written as 

1990 :math:`A_0 \otimes \cdots \otimes A_{n-1}` for matrices 

1991 :math:`A_0, \ldots, A_{n-1}` with shapes given in ``dimensions``, then 

1992 this returns :math:`B_0 \otimes \cdots \otimes B_{n-1}` with 

1993 :math:`B_i = \operatorname{tr}(A_i)`, if ``i in subsystems`` (with 

1994 :math:`i = -1` read as :math:`n-1`), and :math:`B_i = A_i`, otherwise. 

1995 

1996 :param subsystems: A collection of or a single subystem number, indexed 

1997 from zero, corresponding to subsystems that shall be traced over. 

1998 The value :math:`-1` refers to the last subsystem. 

1999 :type subsystems: int or tuple or list 

2000 

2001 :param dimensions: Either an integer :math:`d` so that the subsystems 

2002 are assumed to be all of shape :math:`d \times d`, or a sequence of 

2003 subsystem shapes where an integer :math:`d` within the sequence is 

2004 read as :math:`d \times d`. In any case, the elementwise product 

2005 over all subsystem shapes must equal the expression's shape. 

2006 :type dimensions: int or tuple or list 

2007 

2008 :raises TypeError: If the subsystems do not match the expression or if 

2009 a non-square subsystem is to be traced over. 

2010 :raises IndexError: If the subsystem selection is invalid in any other 

2011 way. 

2012 

2013 :Example: 

2014 

2015 >>> from picos import Constant 

2016 >>> A = Constant("A", range(16), (4, 4)) 

2017 >>> print(A) #doctest: +NORMALIZE_WHITESPACE 

2018 [ 0.00e+00 4.00e+00 8.00e+00 1.20e+01] 

2019 [ 1.00e+00 5.00e+00 9.00e+00 1.30e+01] 

2020 [ 2.00e+00 6.00e+00 1.00e+01 1.40e+01] 

2021 [ 3.00e+00 7.00e+00 1.10e+01 1.50e+01] 

2022 >>> A0 = A.partial_trace(0); A0 

2023 <2×2 Real Constant: A.{tr([2×2])⊗[2×2]}> 

2024 >>> print(A0) #doctest: +NORMALIZE_WHITESPACE 

2025 [ 1.00e+01 1.80e+01] 

2026 [ 1.20e+01 2.00e+01] 

2027 >>> A1 = A.partial_trace(1); A1 

2028 <2×2 Real Constant: A.{[2×2]⊗tr([2×2])}> 

2029 >>> print(A1) #doctest: +NORMALIZE_WHITESPACE 

2030 [ 5.00e+00 2.10e+01] 

2031 [ 9.00e+00 2.50e+01] 

2032 """ 

2033 # Shape of the original matrix. 

2034 m, n = self._shape 

2035 

2036 if isinstance(dimensions, int): 

2037 dimensions = self._square_equal_subsystem_dims(dimensions) 

2038 else: 

2039 dimensions = [ 

2040 (d, d) if isinstance(d, int) else d for d in dimensions] 

2041 

2042 if reduce( 

2043 lambda x, y: (x[0]*y[0], x[1]*y[1]), dimensions) != self._shape: 

2044 raise TypeError("Subsystem dimensions do not match expression.") 

2045 

2046 if isinstance(subsystems, int): 

2047 subsystems = (subsystems,) 

2048 elif not subsystems: 

2049 return self 

2050 

2051 numSys = len(dimensions) 

2052 subsystems = set(numSys - 1 if sys == -1 else sys for sys in subsystems) 

2053 

2054 for sys in subsystems: 

2055 if not isinstance(sys, int): 

2056 raise IndexError("Subsystem indices must be integer, not {}." 

2057 .format(type(sys).__name__)) 

2058 elif sys < 0: 

2059 raise IndexError("Subsystem indices must be nonnegative.") 

2060 elif sys >= numSys: 

2061 raise IndexError( 

2062 "Subsystem index {} out of range for {} systems total." 

2063 .format(sys, numSys)) 

2064 elif dimensions[sys][0] != dimensions[sys][1]: 

2065 raise TypeError( 

2066 "Subsystem index {} refers to a non-square subsystem that " 

2067 "cannot be traced over.".format(sys)) 

2068 

2069 # If all subsystems are traced over, this is the regular trace. 

2070 if len(subsystems) == numSys: 

2071 return self.tr 

2072 

2073 # Prepare sparse T such that T·vec(A) = vec(partial_trace(A)). 

2074 T = [] 

2075 

2076 # Compute factors of T, one for each subsystem being traced over. 

2077 ibh, ibw = m, n 

2078 sysStrings = None 

2079 for sys in range(numSys): 

2080 # Shape of current system. 

2081 p, q = dimensions[sys] 

2082 sysString = glyphs.matrix(glyphs.shape((p, q))) 

2083 

2084 # Heigh/width of "outer" blocks whose relative position is 

2085 # maintained but that are reduced independently to the size of an 

2086 # inner block if the current system is to be traced over. 

2087 obh, obw = ibh, ibw 

2088 

2089 # Height/width of "inner" blocks that are summed if they are 

2090 # main-diagonal blocks of an outer block. 

2091 ibh, ibw = obh // p, obw // q 

2092 

2093 # Only trace over selected subsystems. 

2094 if sys not in subsystems: 

2095 sysStrings = glyphs.kron(sysStrings, sysString) \ 

2096 if sysStrings else sysString 

2097 continue 

2098 else: 

2099 sysStrings = glyphs.kron(sysStrings, glyphs.trace(sysString)) \ 

2100 if sysStrings else glyphs.trace(sysString) 

2101 

2102 # Shape of new matrix. 

2103 assert p == q 

2104 mN, nN = m // p, n // p 

2105 

2106 # Prepare one factor of T. 

2107 oldLen = m * n 

2108 newLen = mN * nN 

2109 V, I, J = [1]*(newLen*p), [], [] 

2110 shape = (newLen, oldLen) 

2111 

2112 # Determine the summands that make up each entry of the new matrix. 

2113 for viN in range(newLen): 

2114 # A column/row pair that represents a scalar in the new matrix 

2115 # and a number of p scalars within different on-diagonal inner 

2116 # blocks but within the same outer block in the old matrix. 

2117 cN, rN = divmod(viN, mN) 

2118 

2119 # Index pair (obi, obj) for outer block in question, row/column 

2120 # pair (ibr, ibc) identifies the scalar within each inner block. 

2121 obi, ibr = divmod(rN, ibh) 

2122 obj, ibc = divmod(cN, ibw) 

2123 

2124 # Collect summands for the current entry of the new matrix; one 

2125 # scalar per on-diagonal inner block. 

2126 for k in range(p): 

2127 rO = obi*obh + k*ibh + ibr 

2128 cO = obj*obw + k*ibw + ibc 

2129 I.append(viN) 

2130 J.append(cO*m + rO) 

2131 

2132 # Store the operator that performs the current subsystem trace. 

2133 T.insert(0, cvxopt.spmatrix(V, I, J, shape, self._typecode)) 

2134 

2135 # Make next iteration work on the new matrix. 

2136 m, n = mN, nN 

2137 

2138 # Multiply out the linear partial trace operator T. 

2139 # TODO: Fast matrix multiplication dynamic program useful here? 

2140 T = reduce(lambda A, B: A*B, T) 

2141 

2142 string = glyphs.ptrace_(self.string, sysStrings) 

2143 shape = (m, n) 

2144 coefs = {mtbs: T * coef for mtbs, coef in self._coefs.items()} 

2145 

2146 return self._basetype(string, shape, coefs) 

2147 

2148 @cached_property 

2149 def tr0(self): 

2150 r"""Expression with the first :math:`2 \times 2` subsystem traced out. 

2151 

2152 Only available for a :math:`2^k \times 2^k` matrix with all subsystems 

2153 of shape :math:`2 \times 2`. Use :meth:`partial_trace` otherwise. 

2154 """ 

2155 return self.partial_trace(subsystems=0) 

2156 

2157 @cached_property 

2158 def tr1(self): 

2159 r"""Expression with the second :math:`2 \times 2` subsystem traced out. 

2160 

2161 Only available for a :math:`2^k \times 2^k` matrix with all subsystems 

2162 of shape :math:`2 \times 2`. Use :meth:`partial_trace` otherwise. 

2163 """ 

2164 return self.partial_trace(subsystems=1) 

2165 

2166 @cached_property 

2167 def tr2(self): 

2168 r"""Expression with the third :math:`2 \times 2` subsystem traced out. 

2169 

2170 Only available for a :math:`2^k \times 2^k` matrix with all subsystems 

2171 of shape :math:`2 \times 2`. Use :meth:`partial_trace` otherwise. 

2172 """ 

2173 return self.partial_trace(subsystems=2) 

2174 

2175 @cached_property 

2176 def tr3(self): 

2177 r"""Expression with the fourth :math:`2 \times 2` subsystem traced out. 

2178 

2179 Only available for a :math:`2^k \times 2^k` matrix with all subsystems 

2180 of shape :math:`2 \times 2`. Use :meth:`partial_trace` otherwise. 

2181 """ 

2182 return self.partial_trace(subsystems=3) 

2183 

2184 @cached_property 

2185 def trl(self): 

2186 r"""Expression with the last :math:`2 \times 2` subsystem traced out. 

2187 

2188 Only available for a :math:`2^k \times 2^k` matrix with all subsystems 

2189 of shape :math:`2 \times 2`. Use :meth:`partial_trace` otherwise. 

2190 """ 

2191 return self.partial_trace(subsystems=-1) 

2192 

2193 @cached_property 

2194 def vec(self): 

2195 """Column-major vectorization of the expression as a column vector. 

2196 

2197 .. note:: 

2198 Given an expression ``A``, ``A.vec`` and ``A[:]`` produce the same 

2199 result (up to its string description) but ``A.vec`` is faster and 

2200 its result is cached. 

2201 

2202 :Example: 

2203 

2204 >>> from picos import Constant 

2205 >>> A = Constant("A", [[1, 2], [3, 4]]) 

2206 >>> A.vec.equals(A[:]) 

2207 True 

2208 >>> A[:] is A[:] 

2209 False 

2210 >>> A.vec is A.vec 

2211 True 

2212 """ 

2213 if self._shape[1] == 1: 

2214 return self 

2215 else: 

2216 return self._basetype( 

2217 glyphs.vec(self.string), (len(self), 1), self._coefs) 

2218 

2219 def dupvec(self, n): 

2220 """Return a (repeated) column-major vectorization of the expression. 

2221 

2222 :param int n: Number of times to duplicate the vectorization. 

2223 

2224 :returns: A column vector. 

2225 

2226 :Example: 

2227 

2228 >>> from picos import Constant 

2229 >>> A = Constant("A", [[1, 2], [3, 4]]) 

2230 >>> A.dupvec(1) is A.vec 

2231 True 

2232 >>> A.dupvec(3).equals(A.vec // A.vec // A.vec) 

2233 True 

2234 """ 

2235 if not isinstance(n, int): 

2236 raise TypeError("Number of copies must be integer.") 

2237 

2238 if n < 1: 

2239 raise ValueError("Number of copies must be positive.") 

2240 

2241 if n == 1: 

2242 return self.vec 

2243 else: 

2244 string = glyphs.vec(glyphs.comma(self.string, n)) 

2245 shape = (len(self)*n, 1) 

2246 coefs = {mtbs: cvxopt_vcat([coef]*n) 

2247 for mtbs, coef in self._coefs.items()} 

2248 

2249 return self._basetype(string, shape, coefs) 

2250 

2251 @cached_property 

2252 def trilvec(self): 

2253 r"""Column-major vectorization of the lower triangular part. 

2254 

2255 :returns: 

2256 A column vector of all elements :math:`A_{ij}` that satisfy 

2257 :math:`i \geq j`. 

2258 

2259 .. note:: 

2260 

2261 If you want a row-major vectorization instead, write ``A.T.triuvec`` 

2262 instead of ``A.trilvec``. 

2263 

2264 :Example: 

2265 

2266 >>> from picos import Constant 

2267 >>> A = Constant("A", [[1, 2], [3, 4], [5, 6]]) 

2268 >>> print(A) 

2269 [ 1.00e+00 2.00e+00] 

2270 [ 3.00e+00 4.00e+00] 

2271 [ 5.00e+00 6.00e+00] 

2272 >>> print(A.trilvec) 

2273 [ 1.00e+00] 

2274 [ 3.00e+00] 

2275 [ 5.00e+00] 

2276 [ 4.00e+00] 

2277 [ 6.00e+00] 

2278 """ 

2279 m, n = self._shape 

2280 

2281 if n == 1: # Column vector or scalar. 

2282 return self 

2283 elif m == 1: # Row vector. 

2284 return self[0] 

2285 

2286 # Build a transformation D such that D·vec(A) = trilvec(A). 

2287 rows = [j*m + i for j in range(n) for i in range(m) if i >= j] 

2288 d = len(rows) 

2289 D = cvxopt.spmatrix([1]*d, range(d), rows, (d, len(self))) 

2290 

2291 string = glyphs.trilvec(self.string) 

2292 shape = (d, 1) 

2293 coefs = {mtbs: D*coef for mtbs, coef in self._coefs.items()} 

2294 

2295 return self._basetype(string, shape, coefs) 

2296 

2297 @cached_property 

2298 def triuvec(self): 

2299 r"""Column-major vectorization of the upper triangular part. 

2300 

2301 :returns: 

2302 A column vector of all elements :math:`A_{ij}` that satisfy 

2303 :math:`i \leq j`. 

2304 

2305 .. note:: 

2306 

2307 If you want a row-major vectorization instead, write ``A.T.trilvec`` 

2308 instead of ``A.triuvec``. 

2309 

2310 :Example: 

2311 

2312 >>> from picos import Constant 

2313 >>> A = Constant("A", [[1, 2, 3], [4, 5, 6]]) 

2314 >>> print(A) 

2315 [ 1.00e+00 2.00e+00 3.00e+00] 

2316 [ 4.00e+00 5.00e+00 6.00e+00] 

2317 >>> print(A.triuvec) 

2318 [ 1.00e+00] 

2319 [ 2.00e+00] 

2320 [ 5.00e+00] 

2321 [ 3.00e+00] 

2322 [ 6.00e+00] 

2323 """ 

2324 m, n = self._shape 

2325 

2326 if m == 1: # Row vector or scalar. 

2327 return self 

2328 elif n == 1: # Column vector. 

2329 return self[0] 

2330 

2331 # Build a transformation D such that D·vec(A) = triuvec(A). 

2332 rows = [j*m + i for j in range(n) for i in range(m) if i <= j] 

2333 d = len(rows) 

2334 D = cvxopt.spmatrix([1]*d, range(d), rows, (d, len(self))) 

2335 

2336 string = glyphs.triuvec(self.string) 

2337 shape = (d, 1) 

2338 coefs = {mtbs: D*coef for mtbs, coef in self._coefs.items()} 

2339 

2340 return self._basetype(string, shape, coefs) 

2341 

2342 @cached_property 

2343 def svec(self): 

2344 """An isometric vectorization of a symmetric or hermitian expression. 

2345 

2346 In the real symmetric case 

2347 

2348 - the vectorization format is precisely the one define in [svec]_, 

2349 - the vectorization is isometric and isomorphic, and 

2350 - this is the same vectorization as used internally by the 

2351 :class:`~.variables.SymmetricVariable` class. 

2352 

2353 In the complex hermitian case 

2354 

2355 - the same format is used, now resulting in a complex vector, 

2356 - the vectorization is isometric but **not** isomorphic as there are 

2357 guaranteed zeros in the imaginary part of the vector, and 

2358 - this is **not** the same vectorization as the isomorphic, 

2359 real-valued one used by :class:`~.variables.HermitianVariable`. 

2360 

2361 The reverse operation is denoted by :attr:`desvec` in either case. 

2362 

2363 :raises TypeError: 

2364 If the expression is not square. 

2365 

2366 :raises ValueError: 

2367 If the expression is not hermitian. 

2368 """ 

2369 if not self.square: 

2370 raise TypeError("Can only compute svec for a square matrix, not for" 

2371 " the {} expression {}.".format(self._shape, self.string)) 

2372 elif not self.hermitian: 

2373 raise ValueError("Cannot compute svec for the non-hermitian " 

2374 "expression {}.".format(self.string)) 

2375 

2376 vec = SymmetricVectorization(self._shape) 

2377 V = vec._full2special 

2378 

2379 string = glyphs.svec(self.string) 

2380 

2381 if self.isreal: 

2382 vec = SymmetricVectorization(self._shape) 

2383 V = vec._full2special 

2384 coefs = {mtbs: V*coef for mtbs, coef in self._coefs.items()} 

2385 result = self._basetype(string, (vec.dim, 1), coefs) 

2386 else: 

2387 # NOTE: We need to reproduce svec using triuvec because, for numeric 

2388 # reasons, SymmetricVectorization averages the elements in the 

2389 # lower and upper triangular part. For symmetric matrices, 

2390 # this is equivalent to the formal definition of svec found in 

2391 # Datorro's book. For hermitian matrices however it is not. 

2392 real_svec = self.real.svec 

2393 imag_svec = 2**0.5 * 1j * self.imag.triuvec 

2394 

2395 result = (real_svec + imag_svec).renamed(string) 

2396 

2397 with unlocked_cached_properties(result): 

2398 result.desvec = self 

2399 

2400 return result 

2401 

2402 @cached_property 

2403 def desvec(self): 

2404 r"""The reverse operation of :attr:`svec`. 

2405 

2406 :raises TypeError: 

2407 If the expression is not a vector or has a dimension outside the 

2408 permissible set 

2409 

2410 .. math:: 

2411 

2412 \left\{ \frac{n(n + 1)}{2} 

2413 \mid n \in \mathbb{Z}_{\geq 1} \right\} 

2414 = \left\{ n \in \mathbb{Z}_{\geq 1} 

2415 \mid \frac{1}{2} \left( \sqrt{8n + 1} - 1 \right) 

2416 \in \mathbb{Z}_{\geq 1} \right\}. 

2417 

2418 :raises ValueError: 

2419 In the case of a complex vector, If the vector is not in the image 

2420 of :attr:`svec`. 

2421 """ 

2422 if 1 not in self.shape: 

2423 raise TypeError( 

2424 "Can only compute desvec for a vector, not for the {} " 

2425 "expression {}.".format(glyphs.shape(self._shape), self.string)) 

2426 

2427 n = 0.5*((8*len(self) + 1)**0.5 - 1) 

2428 if int(n) != n: 

2429 raise TypeError("Cannot compute desvec for the {}-dimensional " 

2430 "vector {} as this size is not a possible outcome of svec." 

2431 .format(len(self), self.string)) 

2432 n = int(n) 

2433 

2434 vec = SymmetricVectorization((n, n)) 

2435 D = vec._special2full 

2436 string = glyphs.desvec(self.string) 

2437 

2438 if self.isreal: 

2439 coefs = {mtbs: D*coef for mtbs, coef in self._coefs.items()} 

2440 result = self._basetype(string, (n, n), coefs) 

2441 

2442 assert result.hermitian 

2443 else: 

2444 # NOTE: While :attr:`svec` performs essentially the same operation 

2445 # for both symmetric and hermitian matrices, we now need to 

2446 # treat the imaginary separately as the imaginary part of a 

2447 # hermitian matrix is skew-symmetric instead of symmetric. 

2448 V, I, J = [], D.I, D.J 

2449 for v, i, j in zip(D.V, I, J): 

2450 V.append(1j*v if i % n <= i // n else -1j*v) 

2451 C = cvxopt.spmatrix(V, I, J, D.size, tc="z") 

2452 

2453 real_desvec = self.real.desvec 

2454 imag_desvec = self._basetype(string, (n, n), { 

2455 mtbs: C*coef for mtbs, coef in self.imag._coefs.items()}) 

2456 

2457 result = (real_desvec + imag_desvec).renamed(string) 

2458 

2459 if not result.hermitian: 

2460 raise ValueError("The complex vector {} is not in the image of " 

2461 "svec. Note that svec is not bijective in the complex case." 

2462 .format(self.string)) 

2463 

2464 with unlocked_cached_properties(result): 

2465 result.svec = self 

2466 

2467 return result 

2468 

2469 def dupdiag(self, n): 

2470 """Return a matrix with the (repeated) expression on the diagonal. 

2471 

2472 Vectorization is performed in column-major order. 

2473 

2474 :param int n: Number of times to duplicate the vectorization. 

2475 """ 

2476 # Vectorize and duplicate the expression. 

2477 vec = self.dupvec(n) 

2478 d = len(vec) 

2479 

2480 # Build a transformation D such that D·vec(A) = vec(diag(vec(A))). 

2481 ones = [1]*d 

2482 D = cvxopt.spdiag(ones)[:] 

2483 D = cvxopt.spmatrix(ones, D.I, range(d), (D.size[0], d)) 

2484 

2485 string = glyphs.diag(vec.string) 

2486 shape = (d, d) 

2487 coefs = {mtbs: D*coef for mtbs, coef in vec._coefs.items()} 

2488 

2489 return self._basetype(string, shape, coefs) 

2490 

2491 @cached_property 

2492 def diag(self): 

2493 """Diagonal matrix with the expression on the main diagonal. 

2494 

2495 Vectorization is performed in column-major order. 

2496 """ 

2497 return self.dupdiag(1) 

2498 

2499 @cached_property 

2500 def maindiag(self): 

2501 """The main diagonal of the expression as a column vector.""" 

2502 if 1 in self._shape: 

2503 return self[0] 

2504 

2505 # Build a transformation D such that D·vec(A) = diag(A). 

2506 step = self._shape[0] + 1 

2507 rows = [i*step for i in range(min(self._shape))] 

2508 d = len(rows) 

2509 D = cvxopt.spmatrix([1]*d, range(d), rows, (d, len(self))) 

2510 

2511 string = glyphs.maindiag(self.string) 

2512 shape = (d, 1) 

2513 coefs = {mtbs: D*coef for mtbs, coef in self._coefs.items()} 

2514 

2515 return self._basetype(string, shape, coefs) 

2516 

2517 def factor_out(self, mutable): 

2518 r"""Factor out a single mutable from a vector expression. 

2519 

2520 If this expression is a column vector :math:`a` that depends on some 

2521 mutable :math:`x` with a trivial internal vectorization format (i.e. 

2522 :class:`~.vectorizations.FullVectorization`), then this method, called 

2523 with :math:`x` as its argument, returns a pair of expressions 

2524 :math:`(a_x, a_0)` such that :math:`a = a_x\operatorname{vec}(x) + a_0`. 

2525 

2526 :returns: 

2527 Two refined :class:`BiaffineExpression` instances that do not depend 

2528 on ``mutable``. 

2529 

2530 :raises TypeError: 

2531 If the expression is not a column vector or if ``mutable`` is not a 

2532 :class:`~.mutable.Mutable` or does not have a trivial vectorization 

2533 format. 

2534 

2535 :raises LookupError: 

2536 If the expression does not depend on ``mutable``. 

2537 

2538 :Example: 

2539 

2540 >>> from picos import RealVariable 

2541 >>> from picos.uncertain import UnitBallPerturbationSet 

2542 >>> x = RealVariable("x", 3) 

2543 >>> z = UnitBallPerturbationSet("z", 3).parameter 

2544 >>> a = ((2*x + 3)^z) + 4*x + 5; a 

2545 <3×1 Uncertain Affine Expression: (2·x + [3])⊙z + 4·x + [5]> 

2546 >>> sorted(a.mutables, key=lambda mtb: mtb.name) 

2547 [<3×1 Real Variable: x>, <3×1 Perturbation: z>] 

2548 >>> az, a0 = a.factor_out(z) 

2549 >>> az 

2550 <3×3 Real Affine Expression: ((2·x + [3])⊙z + 4·x + [5])_z> 

2551 >>> a0 

2552 <3×1 Real Affine Expression: ((2·x + [3])⊙z + 4·x + [5])_0> 

2553 >>> sorted(az.mutables.union(a0.mutables), key=lambda mtb: mtb.name) 

2554 [<3×1 Real Variable: x>] 

2555 >>> (az*z + a0).equals(a) 

2556 True 

2557 """ 

2558 if self._shape[1] != 1: 

2559 raise TypeError( 

2560 "Can only factor out mutables from column vectors but {} has " 

2561 "a shape of {}.".format(self.string, glyphs.shape(self._shape))) 

2562 

2563 mtb = mutable 

2564 

2565 if not isinstance(mtb, Mutable): 

2566 raise TypeError("Can only factor out mutables, not instances of {}." 

2567 .format(type(mtb).__name__)) 

2568 

2569 if not isinstance(mtb._vec, FullVectorization): 

2570 raise TypeError("Can only factor out mutables with a trivial " 

2571 "vectorization format but {} uses {}." 

2572 .format(mtb.name, type(mtb._vec).__name__)) 

2573 

2574 if mtb not in self.mutables: 

2575 raise LookupError("Cannot factor out {} from {} as the latter does " 

2576 "not depend on the former.".format(mtb.name, self.string)) 

2577 

2578 ax_string = glyphs.index(self.string, mtb.name) 

2579 ax_shape = (self._shape[0], mtb.dim) 

2580 

2581 # Recipe in "Robust conic optimization in Python" (Stahlberg 2020). 

2582 ax_coefs = {} 

2583 for mtbs, coef in self._coefs.items(): 

2584 if mtb not in set(mtbs): 

2585 continue 

2586 elif len(mtbs) == 1: 

2587 assert mtbs[0] is mtb 

2588 self._update_coefs(ax_coefs, (), coef[:]) 

2589 else: 

2590 if mtbs[0] is mtb: 

2591 other_mtb = mtbs[1] 

2592 C = coef*cvxopt_K(other_mtb.dim, mtb.dim) 

2593 else: 

2594 assert mtbs[1] is mtb 

2595 other_mtb = mtbs[0] 

2596 C = coef 

2597 

2598 d = other_mtb.dim 

2599 C = cvxopt.sparse([C[:, i*d:(i+1)*d] for i in range(mtb.dim)]) 

2600 self._update_coefs(ax_coefs, (other_mtb,), C) 

2601 

2602 ax = self._basetype(ax_string, ax_shape, ax_coefs) 

2603 

2604 a0_string = glyphs.index(self.string, 0) 

2605 a0_shape = self._shape 

2606 a0_coefs = {M: C for M, C in self._coefs.items() if mtb not in set(M)} 

2607 a0 = self._basetype(a0_string, a0_shape, a0_coefs) 

2608 

2609 return ax.refined, a0.refined 

2610 

2611 # -------------------------------------------------------------------------- 

2612 # Backwards compatibility methods. 

2613 # -------------------------------------------------------------------------- 

2614 

2615 @classmethod 

2616 @deprecated("2.0", useInstead="from_constant") 

2617 def fromScalar(cls, scalar): 

2618 """Create a class instance from a numeric scalar.""" 

2619 return cls.from_constant(scalar, (1, 1)) 

2620 

2621 @classmethod 

2622 @deprecated("2.0", useInstead="from_constant") 

2623 def fromMatrix(cls, matrix, size=None): 

2624 """Create a class instance from a numeric matrix.""" 

2625 return cls.from_constant(matrix, size) 

2626 

2627 @deprecated("2.0", useInstead="object.__xor__") 

2628 def hadamard(self, fact): 

2629 """Denote the elementwise (or Hadamard) product.""" 

2630 return self.__xor__(fact) 

2631 

2632 @deprecated("2.0", useInstead="~.expression.Expression.constant") 

2633 def isconstant(self): 

2634 """Whether the expression involves no mutables.""" 

2635 return self.constant 

2636 

2637 @deprecated("2.0", useInstead="equals") 

2638 def same_as(self, other): 

2639 """Check mathematical equality with another affine expression.""" 

2640 return self.equals(other) 

2641 

2642 @deprecated("2.0", useInstead="T") 

2643 def transpose(self): 

2644 """Return the matrix transpose.""" 

2645 return self.T 

2646 

2647 @cached_property 

2648 @deprecated("2.0", useInstead="partial_transpose", decoratorLevel=1) 

2649 def Tx(self): 

2650 """Auto-detect few subsystems of same shape and transpose the last.""" 

2651 m, n = self._shape 

2652 dims = None 

2653 

2654 for k in range(2, int(math.log(min(m, n), 2)) + 1): 

2655 p, q = int(round(m**(1.0/k))), int(round(n**(1.0/k))) 

2656 if m == p**k and n == q**k: 

2657 dims = ((p, q),)*k 

2658 break 

2659 

2660 if dims: 

2661 return self.partial_transpose(subsystems=-1, dimensions=dims) 

2662 else: 

2663 raise RuntimeError("Failed to auto-detect subsystem dimensions for " 

2664 "partial transposition: Only supported for {} matrices, {}." 

2665 .format(glyphs.shape( 

2666 (glyphs.power("m", "k"), glyphs.power("n", "k"))), 

2667 glyphs.ge("k", 2))) 

2668 

2669 @deprecated("2.0", useInstead="conj") 

2670 def conjugate(self): 

2671 """Return the complex conjugate.""" 

2672 return self.conj 

2673 

2674 @deprecated("2.0", useInstead="H") 

2675 def Htranspose(self): 

2676 """Return the conjugate (or Hermitian) transpose.""" 

2677 return self.H 

2678 

2679 @deprecated("2.0", reason="PICOS expressions are now immutable.") 

2680 def copy(self): 

2681 """Return a deep copy of the expression.""" 

2682 from copy import copy as cp 

2683 

2684 return self._basetype(cp(self._symbStr), self._shape, 

2685 {mtbs: cp(coef) for mtbs, coef in self._coefs.items()}) 

2686 

2687 @deprecated("2.0", reason="PICOS expressions are now immutable.") 

2688 def soft_copy(self): 

2689 """Return a shallow copy of the expression.""" 

2690 return self._basetype(self._symbStr, self._shape, self._coefs) 

2691 

2692 

2693# -------------------------------------- 

2694__all__ = api_end(_API_START, globals())