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

1373 """Return the expression reshaped in column-major order. 

1374 

1375 :Example: 

1376 

1377 >>> from picos import Constant 

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

1379 >>> print(C) 

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

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

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

1383 [ 0.00e+00 3.00e+00] 

1384 [ 1.00e+00 4.00e+00] 

1385 [ 2.00e+00 5.00e+00] 

1386 """ 

1387 shape = load_shape(shape, wildcards=True) 

1388 

1389 if shape == self._shape: 

1390 return self 

1391 elif shape == (None, None): 

1392 return self 

1393 

1394 length = len(self) 

1395 

1396 if shape[0] is None: 

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

1398 elif shape[1] is None: 

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

1400 

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

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

1403 

1404 string = glyphs.reshaped(self.string, glyphs.shape(shape)) 

1405 return self._basetype(string, shape, self._coefs) 

1406 

1407 def broadcasted(self, shape): 

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

1409 

1410 :Example: 

1411 

1412 >>> from picos import Constant 

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

1414 >>> print(C) 

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

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

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

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

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

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

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

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

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

1424 """ 

1425 shape = load_shape(shape, wildcards=True) 

1426 shape = blend_shapes(shape, self._shape) 

1427 

1428 if shape == self._shape: 

1429 return self 

1430 

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

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

1433 

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

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

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

1437 

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

1439 string = glyphs.matrix(self.string) 

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

1441 

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

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

1444 

1445 def reshaped_or_broadcasted(self, shape): 

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

1447 

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

1449 may not contain a wildcard character. 

1450 

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

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

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

1454 """ 

1455 shape = load_shape(shape, wildcards=False) 

1456 

1457 try: 

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

1459 return self.reshaped(shape) 

1460 else: 

1461 return self.broadcasted(shape) 

1462 except ValueError: 

1463 raise ValueError( 

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

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

1466 

1467 @cached_property 

1468 def hermitianized(self): 

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

1470 

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

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

1473 

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

1475 subspace of symmetric matrices. 

1476 """ 

1477 if not self.square: 

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

1479 

1480 return (self + self.H)/2 

1481 

1482 @cached_property 

1483 def real(self): 

1484 """Real part of the expression.""" 

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

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

1487 

1488 @cached_property 

1489 def imag(self): 

1490 """Imaginary part of the expression.""" 

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

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

1493 

1494 @cached_property 

1495 def bilin(self): 

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

1497 return self._basetype( 

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

1499 

1500 @cached_property 

1501 def lin(self): 

1502 """Linear part of the expression.""" 

1503 return self._basetype( 

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

1505 

1506 @cached_property 

1507 def noncst(self): 

1508 """Nonconstant part of the expression.""" 

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

1510 

1511 return self._basetype( 

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

1513 

1514 @cached_property 

1515 def cst(self): 

1516 """Constant part of the expression.""" 

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

1518 

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

1520 

1521 @cached_selfinverse_property 

1522 def T(self): 

1523 """Matrix transpose.""" 

1524 if len(self) == 1: 

1525 return self 

1526 

1527 m, n = self._shape 

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

1529 

1530 string = glyphs.transp(self.string) 

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

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

1533 

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

1535 

1536 @cached_selfinverse_property 

1537 def conj(self): 

1538 """Complex conjugate.""" 

1539 string = glyphs.conj(self.string) 

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

1541 

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

1543 

1544 @cached_selfinverse_property 

1545 def H(self): 

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

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

1548 

1549 def _square_equal_subsystem_dims(self, diagLen): 

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

1551 m, n = self._shape 

1552 k = math.log(m, diagLen) 

1553 

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

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

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

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

1558 

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

1560 

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

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

1563 

1564 If the expression can be written as 

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

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

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

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

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

1570 

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

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

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

1574 :type subsystems: int or tuple or list 

1575 

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

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

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

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

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

1581 :type dimensions: int or tuple or list 

1582 

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

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

1585 

1586 :Example: 

1587 

1588 >>> from picos import Constant 

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

1609 """ 

1610 m, n = self._shape 

1611 

1612 if isinstance(dimensions, int): 

1613 dimensions = self._square_equal_subsystem_dims(dimensions) 

1614 else: 

1615 dimensions = [ 

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

1617 

1618 if reduce( 

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

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

1621 

1622 if isinstance(subsystems, int): 

1623 subsystems = (subsystems,) 

1624 elif not subsystems: 

1625 return self 

1626 

1627 numSys = len(dimensions) 

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

1629 

1630 for sys in subsystems: 

1631 if not isinstance(sys, int): 

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

1633 .format(type(sys).__name__)) 

1634 elif sys < 0: 

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

1636 elif sys >= numSys: 

1637 raise IndexError( 

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

1639 .format(sys, numSys)) 

1640 

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

1642 if len(subsystems) == numSys: 

1643 return self.T 

1644 

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

1646 d = m * n 

1647 V = [1]*d 

1648 I = range(d) 

1649 J = cvxopt.matrix(I) 

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

1651 obh, obw = 1, 1 

1652 sysStrings = None 

1653 

1654 # Apply transpositions starting with the rightmost Kronecker factor. 

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

1656 # Shape of current system. 

1657 p, q = dimensions[sys] 

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

1659 

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

1661 ibh, ibw = obh, obw 

1662 

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

1664 # maintained but that are subject to transposition independently. 

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

1666 obh *= p 

1667 obw *= q 

1668 

1669 # Only transpose selected subsystems. 

1670 if sys not in subsystems: 

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

1672 if sysStrings else sysString 

1673 continue 

1674 else: 

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

1676 if sysStrings else glyphs.transp(sysString) 

1677 

1678 # Shape of outer blocks after transposition. 

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

1680 

1681 # Shape of full matrix after transposition. 

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

1683 

1684 for vi in I: 

1685 # Full matrix column and row. 

1686 c, r = divmod(vi, m) 

1687 

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

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

1690 obi, obr = divmod(r, obh) 

1691 obj, obc = divmod(c, obw) 

1692 

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

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

1695 ibi, ibr = divmod(obr, ibh) 

1696 ibj, ibc = divmod(obc, ibw) 

1697 

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

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

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

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

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

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

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

1705 

1706 # Prepare the transposition. 

1707 T[tvi] = J[vi] 

1708 

1709 # Apply the transposition. 

1710 J, T = T, J 

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

1712 

1713 # Finalize the partial commutation matrix. 

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

1715 

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

1717 shape = (m, n) 

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

1719 

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

1721 

1722 @cached_property 

1723 def T0(self): 

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

1725 

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

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

1728 """ 

1729 return self.partial_transpose(subsystems=0) 

1730 

1731 @cached_property 

1732 def T1(self): 

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

1734 

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

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

1737 """ 

1738 return self.partial_transpose(subsystems=1) 

1739 

1740 @cached_property 

1741 def T2(self): 

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

1743 

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

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

1746 """ 

1747 return self.partial_transpose(subsystems=2) 

1748 

1749 @cached_property 

1750 def T3(self): 

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

1752 

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

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

1755 """ 

1756 return self.partial_transpose(subsystems=3) 

1757 

1758 @cached_property 

1759 def Tl(self): 

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

1761 

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

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

1764 """ 

1765 return self.partial_transpose(subsystems=-1) 

1766 

1767 @staticmethod 

1768 def _reindex_F(indices, source, destination): 

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

1770 new = [] 

1771 offset = 0 

1772 factor = 1 

1773 

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

1775 offset += factor*index 

1776 factor *= base 

1777 

1778 for base in destination: 

1779 offset, remainder = divmod(offset, base) 

1780 new.append(remainder) 

1781 

1782 return tuple(new) 

1783 

1784 @staticmethod 

1785 def _reindex_C(indices, source, destination): 

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

1787 new = [] 

1788 offset = 0 

1789 factor = 1 

1790 

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

1792 offset += factor*index 

1793 factor *= base 

1794 

1795 for base in reversed(destination): 

1796 offset, remainder = divmod(offset, base) 

1797 new.insert(0, remainder) 

1798 

1799 return tuple(new) 

1800 

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

1802 """Return the reshuffled or realigned expression. 

1803 

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

1805 the following sequence of operations: 

1806 

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

1808 according to ``order``. 

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

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

1811 according to ``order``. 

1812 

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

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

1815 

1816 .. code:: 

1817 

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

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

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

1821 reshuffled = numpy.einsum(P, reshuffled) 

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

1823 

1824 :param permutation: 

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

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

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

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

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

1830 should be swapped. 

1831 :type permutation: 

1832 str or tuple or list 

1833 

1834 :param dimensions: 

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

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

1837 hypercubic and the number of dimensions is inferred from the 

1838 ``permutation`` argument. 

1839 :type dimensions: 

1840 None or tuple or list 

1841 

1842 :param str order: 

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

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

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

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

1847 

1848 :Example: 

1849 

1850 >>> from picos import Constant 

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

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

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

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

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

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

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

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

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

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

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

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

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

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

1865 True 

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

1867 True 

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

1869 True 

1870 """ 

1871 m, n = self._shape 

1872 mn = m*n 

1873 

1874 # Load the permutation. 

1875 ordered = sorted(permutation) 

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

1877 

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

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

1880 "contain duplicate elements.") 

1881 

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

1883 

1884 numDims = len(P) 

1885 

1886 # Load the dimensions. 

1887 guessDimensions = dimensions is None 

1888 

1889 if guessDimensions: 

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

1891 else: 

1892 if len(dimensions) != numDims: 

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

1894 "number of dimensions.") 

1895 

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

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

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

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

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

1901 

1902 # Load the indexing order. 

1903 if order not in "FC": 

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

1905 

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

1907 

1908 # Nothing to do for the neutral permutation. 

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

1910 return self 

1911 

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

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

1914 for i in range(mn): 

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

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

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

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

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

1920 I.append(k*m + j) 

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

1922 

1923 # Create the string. 

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

1925 

1926 if not guessDimensions: 

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

1928 

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

1930 

1931 # Finalize the new expression. 

1932 shape = (m, n) 

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

1934 

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

1936 

1937 @cached_property 

1938 def sum(self): 

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

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

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

1942 return (self | 1) 

1943 

1944 @cached_property 

1945 def tr(self): 

1946 """Trace of a square expression.""" 

1947 if not self.square: 

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

1949 .format(self.string)) 

1950 

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

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

1953 return (self | "I") 

1954 

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

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

1957 

1958 If the expression can be written as 

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

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

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

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

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

1964 

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

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

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

1968 :type subsystems: int or tuple or list 

1969 

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

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

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

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

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

1975 :type dimensions: int or tuple or list 

1976 

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

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

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

1980 way. 

1981 

1982 :Example: 

1983 

1984 >>> from picos import Constant 

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

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

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

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

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

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

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

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

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

1994 [ 1.00e+01 1.80e+01] 

1995 [ 1.20e+01 2.00e+01] 

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

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

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

1999 [ 5.00e+00 2.10e+01] 

2000 [ 9.00e+00 2.50e+01] 

2001 """ 

2002 # Shape of the original matrix. 

2003 m, n = self._shape 

2004 

2005 if isinstance(dimensions, int): 

2006 dimensions = self._square_equal_subsystem_dims(dimensions) 

2007 else: 

2008 dimensions = [ 

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

2010 

2011 if reduce( 

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

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

2014 

2015 if isinstance(subsystems, int): 

2016 subsystems = (subsystems,) 

2017 elif not subsystems: 

2018 return self 

2019 

2020 numSys = len(dimensions) 

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

2022 

2023 for sys in subsystems: 

2024 if not isinstance(sys, int): 

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

2026 .format(type(sys).__name__)) 

2027 elif sys < 0: 

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

2029 elif sys >= numSys: 

2030 raise IndexError( 

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

2032 .format(sys, numSys)) 

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

2034 raise TypeError( 

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

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

2037 

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

2039 if len(subsystems) == numSys: 

2040 return self.tr 

2041 

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

2043 T = [] 

2044 

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

2046 ibh, ibw = m, n 

2047 sysStrings = None 

2048 for sys in range(numSys): 

2049 # Shape of current system. 

2050 p, q = dimensions[sys] 

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

2052 

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

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

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

2056 obh, obw = ibh, ibw 

2057 

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

2059 # main-diagonal blocks of an outer block. 

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

2061 

2062 # Only trace over selected subsystems. 

2063 if sys not in subsystems: 

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

2065 if sysStrings else sysString 

2066 continue 

2067 else: 

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

2069 if sysStrings else glyphs.trace(sysString) 

2070 

2071 # Shape of new matrix. 

2072 assert p == q 

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

2074 

2075 # Prepare one factor of T. 

2076 oldLen = m * n 

2077 newLen = mN * nN 

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

2079 shape = (newLen, oldLen) 

2080 

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

2082 for viN in range(newLen): 

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

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

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

2086 cN, rN = divmod(viN, mN) 

2087 

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

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

2090 obi, ibr = divmod(rN, ibh) 

2091 obj, ibc = divmod(cN, ibw) 

2092 

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

2094 # scalar per on-diagonal inner block. 

2095 for k in range(p): 

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

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

2098 I.append(viN) 

2099 J.append(cO*m + rO) 

2100 

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

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

2103 

2104 # Make next iteration work on the new matrix. 

2105 m, n = mN, nN 

2106 

2107 # Multiply out the linear partial trace operator T. 

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

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

2110 

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

2112 shape = (m, n) 

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

2114 

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

2116 

2117 @cached_property 

2118 def tr0(self): 

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

2120 

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

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

2123 """ 

2124 return self.partial_trace(subsystems=0) 

2125 

2126 @cached_property 

2127 def tr1(self): 

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

2129 

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

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

2132 """ 

2133 return self.partial_trace(subsystems=1) 

2134 

2135 @cached_property 

2136 def tr2(self): 

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

2138 

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

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

2141 """ 

2142 return self.partial_trace(subsystems=2) 

2143 

2144 @cached_property 

2145 def tr3(self): 

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

2147 

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

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

2150 """ 

2151 return self.partial_trace(subsystems=3) 

2152 

2153 @cached_property 

2154 def trl(self): 

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

2156 

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

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

2159 """ 

2160 return self.partial_trace(subsystems=-1) 

2161 

2162 @cached_property 

2163 def vec(self): 

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

2165 

2166 .. note:: 

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

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

2169 its result is cached. 

2170 

2171 :Example: 

2172 

2173 >>> from picos import Constant 

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

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

2176 True 

2177 >>> A[:] is A[:] 

2178 False 

2179 >>> A.vec is A.vec 

2180 True 

2181 """ 

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

2183 return self 

2184 else: 

2185 return self._basetype( 

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

2187 

2188 def dupvec(self, n): 

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

2190 

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

2192 

2193 :returns: A column vector. 

2194 

2195 :Example: 

2196 

2197 >>> from picos import Constant 

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

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

2200 True 

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

2202 True 

2203 """ 

2204 if not isinstance(n, int): 

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

2206 

2207 if n < 1: 

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

2209 

2210 if n == 1: 

2211 return self.vec 

2212 else: 

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

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

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

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

2217 

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

2219 

2220 @cached_property 

2221 def trilvec(self): 

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

2223 

2224 :returns: 

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

2226 :math:`i \geq j`. 

2227 

2228 .. note:: 

2229 

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

2231 instead of ``A.trilvec``. 

2232 

2233 :Example: 

2234 

2235 >>> from picos import Constant 

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

2237 >>> print(A) 

2238 [ 1.00e+00 2.00e+00] 

2239 [ 3.00e+00 4.00e+00] 

2240 [ 5.00e+00 6.00e+00] 

2241 >>> print(A.trilvec) 

2242 [ 1.00e+00] 

2243 [ 3.00e+00] 

2244 [ 5.00e+00] 

2245 [ 4.00e+00] 

2246 [ 6.00e+00] 

2247 """ 

2248 m, n = self._shape 

2249 

2250 if n == 1: # Column vector or scalar. 

2251 return self 

2252 elif m == 1: # Row vector. 

2253 return self[0] 

2254 

2255 # Build a transformation D such that D·vec(A) = trilvec(A). 

2256 rows = [j*m + i for j in range(n) for i in range(m) if i >= j] 

2257 d = len(rows) 

2258 D = cvxopt.spmatrix([1]*d, range(d), rows, (d, len(self))) 

2259 

2260 string = glyphs.trilvec(self.string) 

2261 shape = (d, 1) 

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

2263 

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

2265 

2266 @cached_property 

2267 def triuvec(self): 

2268 r"""Column-major vectorization of the upper triangular part. 

2269 

2270 :returns: 

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

2272 :math:`i \leq j`. 

2273 

2274 .. note:: 

2275 

2276 If you want a row-major vectorization instead, write ``A.T.trilvec`` 

2277 instead of ``A.triuvec``. 

2278 

2279 :Example: 

2280 

2281 >>> from picos import Constant 

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

2283 >>> print(A) 

2284 [ 1.00e+00 2.00e+00 3.00e+00] 

2285 [ 4.00e+00 5.00e+00 6.00e+00] 

2286 >>> print(A.triuvec) 

2287 [ 1.00e+00] 

2288 [ 2.00e+00] 

2289 [ 5.00e+00] 

2290 [ 3.00e+00] 

2291 [ 6.00e+00] 

2292 """ 

2293 m, n = self._shape 

2294 

2295 if m == 1: # Row vector or scalar. 

2296 return self 

2297 elif n == 1: # Column vector. 

2298 return self[0] 

2299 

2300 # Build a transformation D such that D·vec(A) = triuvec(A). 

2301 rows = [j*m + i for j in range(n) for i in range(m) if i <= j] 

2302 d = len(rows) 

2303 D = cvxopt.spmatrix([1]*d, range(d), rows, (d, len(self))) 

2304 

2305 string = glyphs.triuvec(self.string) 

2306 shape = (d, 1) 

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

2308 

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

2310 

2311 @cached_property 

2312 def svec(self): 

2313 """An isometric vectorization of a symmetric or hermitian expression. 

2314 

2315 In the real symmetric case 

2316 

2317 - the vectorization format is precisely the one define in [svec]_, 

2318 - the vectorization is isometric and isomorphic, and 

2319 - this is the same vectorization as used internally by the 

2320 :class:`~.variables.SymmetricVariable` class. 

2321 

2322 In the complex hermitian case 

2323 

2324 - the same format is used, now resulting in a complex vector, 

2325 - the vectorization is isometric but **not** isomorphic as there are 

2326 guaranteed zeros in the imaginary part of the vector, and 

2327 - this is **not** the same vectorization as the isomorphic, 

2328 real-valued one used by :class:`~.variables.HermitianVariable`. 

2329 

2330 The reverse operation is denoted by :attr:`desvec` in either case. 

2331 

2332 :raises TypeError: 

2333 If the expression is not square. 

2334 

2335 :raises ValueError: 

2336 If the expression is not hermitian. 

2337 """ 

2338 if not self.square: 

2339 raise TypeError("Can only compute svec for a square matrix, not for" 

2340 " the {} expression {}.".format(self._shape, self.string)) 

2341 elif not self.hermitian: 

2342 raise ValueError("Cannot compute svec for the non-hermitian " 

2343 "expression {}.".format(self.string)) 

2344 

2345 vec = SymmetricVectorization(self._shape) 

2346 V = vec._full2special 

2347 

2348 string = glyphs.svec(self.string) 

2349 

2350 if self.isreal: 

2351 vec = SymmetricVectorization(self._shape) 

2352 V = vec._full2special 

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

2354 result = self._basetype(string, (vec.dim, 1), coefs) 

2355 else: 

2356 # NOTE: We need to reproduce svec using triuvec because, for numeric 

2357 # reasons, SymmetricVectorization averages the elements in the 

2358 # lower and upper triangular part. For symmetric matrices, 

2359 # this is equivalent to the formal definition of svec found in 

2360 # Datorro's book. For hermitian matrices however it is not. 

2361 real_svec = self.real.svec 

2362 imag_svec = 2**0.5 * 1j * self.imag.triuvec 

2363 

2364 result = (real_svec + imag_svec).renamed(string) 

2365 

2366 with unlocked_cached_properties(result): 

2367 result.desvec = self 

2368 

2369 return result 

2370 

2371 @cached_property 

2372 def desvec(self): 

2373 r"""The reverse operation of :attr:`svec`. 

2374 

2375 :raises TypeError: 

2376 If the expression is not a vector or has a dimension outside the 

2377 permissible set 

2378 

2379 .. math:: 

2380 

2381 \left\{ \frac{n(n + 1)}{2} 

2382 \mid n \in \mathbb{Z}_{\geq 1} \right\} 

2383 = \left\{ n \in \mathbb{Z}_{\geq 1} 

2384 \mid \frac{1}{2} \left( \sqrt{8n + 1} - 1 \right) 

2385 \in \mathbb{Z}_{\geq 1} \right\}. 

2386 

2387 :raises ValueError: 

2388 In the case of a complex vector, If the vector is not in the image 

2389 of :attr:`svec`. 

2390 """ 

2391 if 1 not in self.shape: 

2392 raise TypeError( 

2393 "Can only compute desvec for a vector, not for the {} " 

2394 "expression {}.".format(glyphs.shape(self._shape), self.string)) 

2395 

2396 n = 0.5*((8*len(self) + 1)**0.5 - 1) 

2397 if int(n) != n: 

2398 raise TypeError("Cannot compute desvec for the {}-dimensional " 

2399 "vector {} as this size is not a possible outcome of svec." 

2400 .format(len(self), self.string)) 

2401 n = int(n) 

2402 

2403 vec = SymmetricVectorization((n, n)) 

2404 D = vec._special2full 

2405 string = glyphs.desvec(self.string) 

2406 

2407 if self.isreal: 

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

2409 result = self._basetype(string, (n, n), coefs) 

2410 

2411 assert result.hermitian 

2412 else: 

2413 # NOTE: While :attr:`svec` performs essentially the same operation 

2414 # for both symmetric and hermitian matrices, we now need to 

2415 # treat the imaginary separately as the imaginary part of a 

2416 # hermitian matrix is skew-symmetric instead of symmetric. 

2417 V, I, J = [], D.I, D.J 

2418 for v, i, j in zip(D.V, I, J): 

2419 V.append(1j*v if i % n <= i // n else -1j*v) 

2420 C = cvxopt.spmatrix(V, I, J, D.size, tc="z") 

2421 

2422 real_desvec = self.real.desvec 

2423 imag_desvec = self._basetype(string, (n, n), { 

2424 mtbs: C*coef for mtbs, coef in self.imag._coefs.items()}) 

2425 

2426 result = (real_desvec + imag_desvec).renamed(string) 

2427 

2428 if not result.hermitian: 

2429 raise ValueError("The complex vector {} is not in the image of " 

2430 "svec. Note that svec is not bijective in the complex case." 

2431 .format(self.string)) 

2432 

2433 with unlocked_cached_properties(result): 

2434 result.svec = self 

2435 

2436 return result 

2437 

2438 def dupdiag(self, n): 

2439 """Return a matrix with the (repeated) expression on the diagonal. 

2440 

2441 Vectorization is performed in column-major order. 

2442 

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

2444 """ 

2445 # Vectorize and duplicate the expression. 

2446 vec = self.dupvec(n) 

2447 d = len(vec) 

2448 

2449 # Build a transformation D such that D·vec(A) = vec(diag(vec(A))). 

2450 ones = [1]*d 

2451 D = cvxopt.spdiag(ones)[:] 

2452 D = cvxopt.spmatrix(ones, D.I, range(d), (D.size[0], d)) 

2453 

2454 string = glyphs.diag(vec.string) 

2455 shape = (d, d) 

2456 coefs = {mtbs: D*coef for mtbs, coef in vec._coefs.items()} 

2457 

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

2459 

2460 @cached_property 

2461 def diag(self): 

2462 """Diagonal matrix with the expression on the main diagonal. 

2463 

2464 Vectorization is performed in column-major order. 

2465 """ 

2466 return self.dupdiag(1) 

2467 

2468 @cached_property 

2469 def maindiag(self): 

2470 """The main diagonal of the expression as a column vector.""" 

2471 if 1 in self._shape: 

2472 return self[0] 

2473 

2474 # Build a transformation D such that D·vec(A) = diag(A). 

2475 step = self._shape[0] + 1 

2476 rows = [i*step for i in range(min(self._shape))] 

2477 d = len(rows) 

2478 D = cvxopt.spmatrix([1]*d, range(d), rows, (d, len(self))) 

2479 

2480 string = glyphs.maindiag(self.string) 

2481 shape = (d, 1) 

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

2483 

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

2485 

2486 def factor_out(self, mutable): 

2487 r"""Factor out a single mutable from a vector expression. 

2488 

2489 If this expression is a column vector :math:`a` that depends on some 

2490 mutable :math:`x` with a trivial internal vectorization format (i.e. 

2491 :class:`~.vectorizations.FullVectorization`), then this method, called 

2492 with :math:`x` as its argument, returns a pair of expressions 

2493 :math:`(a_x, a_0)` such that :math:`a = a_x\operatorname{vec}(x) + a_0`. 

2494 

2495 :returns: 

2496 Two refined :class:`BiaffineExpression` instances that do not depend 

2497 on ``mutable``. 

2498 

2499 :raises TypeError: 

2500 If the expression is not a column vector or if ``mutable`` is not a 

2501 :class:`~.mutable.Mutable` or does not have a trivial vectorization 

2502 format. 

2503 

2504 :raises LookupError: 

2505 If the expression does not depend on ``mutable``. 

2506 

2507 :Example: 

2508 

2509 >>> from picos import RealVariable 

2510 >>> from picos.uncertain import UnitBallPerturbationSet 

2511 >>> x = RealVariable("x", 3) 

2512 >>> z = UnitBallPerturbationSet("z", 3).parameter 

2513 >>> a = ((2*x + 3)^z) + 4*x + 5; a 

2514 <3×1 Uncertain Affine Expression: (2·x + [3])⊙z + 4·x + [5]> 

2515 >>> sorted(a.mutables, key=lambda mtb: mtb.name) 

2516 [<3×1 Real Variable: x>, <3×1 Perturbation: z>] 

2517 >>> az, a0 = a.factor_out(z) 

2518 >>> az 

2519 <3×3 Real Affine Expression: ((2·x + [3])⊙z + 4·x + [5])_z> 

2520 >>> a0 

2521 <3×1 Real Affine Expression: ((2·x + [3])⊙z + 4·x + [5])_0> 

2522 >>> sorted(az.mutables.union(a0.mutables), key=lambda mtb: mtb.name) 

2523 [<3×1 Real Variable: x>] 

2524 >>> (az*z + a0).equals(a) 

2525 True 

2526 """ 

2527 if self._shape[1] != 1: 

2528 raise TypeError( 

2529 "Can only factor out mutables from column vectors but {} has " 

2530 "a shape of {}.".format(self.string, glyphs.shape(self._shape))) 

2531 

2532 mtb = mutable 

2533 

2534 if not isinstance(mtb, Mutable): 

2535 raise TypeError("Can only factor out mutables, not instances of {}." 

2536 .format(type(mtb).__name__)) 

2537 

2538 if not isinstance(mtb._vec, FullVectorization): 

2539 raise TypeError("Can only factor out mutables with a trivial " 

2540 "vectorization format but {} uses {}." 

2541 .format(mtb.name, type(mtb._vec).__name__)) 

2542 

2543 if mtb not in self.mutables: 

2544 raise LookupError("Cannot factor out {} from {} as the latter does " 

2545 "not depend on the former.".format(mtb.name, self.string)) 

2546 

2547 ax_string = glyphs.index(self.string, mtb.name) 

2548 ax_shape = (self._shape[0], mtb.dim) 

2549 

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

2551 ax_coefs = {} 

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

2553 if mtb not in set(mtbs): 

2554 continue 

2555 elif len(mtbs) == 1: 

2556 assert mtbs[0] is mtb 

2557 self._update_coefs(ax_coefs, (), coef[:]) 

2558 else: 

2559 if mtbs[0] is mtb: 

2560 other_mtb = mtbs[1] 

2561 C = coef*cvxopt_K(other_mtb.dim, mtb.dim) 

2562 else: 

2563 assert mtbs[1] is mtb 

2564 other_mtb = mtbs[0] 

2565 C = coef 

2566 

2567 d = other_mtb.dim 

2568 C = cvxopt.sparse([C[:, i*d:(i+1)*d] for i in range(mtb.dim)]) 

2569 self._update_coefs(ax_coefs, (other_mtb,), C) 

2570 

2571 ax = self._basetype(ax_string, ax_shape, ax_coefs) 

2572 

2573 a0_string = glyphs.index(self.string, 0) 

2574 a0_shape = self._shape 

2575 a0_coefs = {M: C for M, C in self._coefs.items() if mtb not in set(M)} 

2576 a0 = self._basetype(a0_string, a0_shape, a0_coefs) 

2577 

2578 return ax.refined, a0.refined 

2579 

2580 # -------------------------------------------------------------------------- 

2581 # Backwards compatibility methods. 

2582 # -------------------------------------------------------------------------- 

2583 

2584 @classmethod 

2585 @deprecated("2.0", useInstead="from_constant") 

2586 def fromScalar(cls, scalar): 

2587 """Create a class instance from a numeric scalar.""" 

2588 return cls.from_constant(scalar, (1, 1)) 

2589 

2590 @classmethod 

2591 @deprecated("2.0", useInstead="from_constant") 

2592 def fromMatrix(cls, matrix, size=None): 

2593 """Create a class instance from a numeric matrix.""" 

2594 return cls.from_constant(matrix, size) 

2595 

2596 @deprecated("2.0", useInstead="object.__xor__") 

2597 def hadamard(self, fact): 

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

2599 return self.__xor__(fact) 

2600 

2601 @deprecated("2.0", useInstead="~.expression.Expression.constant") 

2602 def isconstant(self): 

2603 """Whether the expression involves no mutables.""" 

2604 return self.constant 

2605 

2606 @deprecated("2.0", useInstead="equals") 

2607 def same_as(self, other): 

2608 """Check mathematical equality with another affine expression.""" 

2609 return self.equals(other) 

2610 

2611 @deprecated("2.0", useInstead="T") 

2612 def transpose(self): 

2613 """Return the matrix transpose.""" 

2614 return self.T 

2615 

2616 @cached_property 

2617 @deprecated("2.0", useInstead="partial_transpose", decoratorLevel=1) 

2618 def Tx(self): 

2619 """Auto-detect few subsystems of same shape and transpose the last.""" 

2620 m, n = self._shape 

2621 dims = None 

2622 

2623 for k in range(2, int(math.log(min(m, n), 2)) + 1): 

2624 p, q = int(round(m**(1.0/k))), int(round(n**(1.0/k))) 

2625 if m == p**k and n == q**k: 

2626 dims = ((p, q),)*k 

2627 break 

2628 

2629 if dims: 

2630 return self.partial_transpose(subsystems=-1, dimensions=dims) 

2631 else: 

2632 raise RuntimeError("Failed to auto-detect subsystem dimensions for " 

2633 "partial transposition: Only supported for {} matrices, {}." 

2634 .format(glyphs.shape( 

2635 (glyphs.power("m", "k"), glyphs.power("n", "k"))), 

2636 glyphs.ge("k", 2))) 

2637 

2638 @deprecated("2.0", useInstead="conj") 

2639 def conjugate(self): 

2640 """Return the complex conjugate.""" 

2641 return self.conj 

2642 

2643 @deprecated("2.0", useInstead="H") 

2644 def Htranspose(self): 

2645 """Return the conjugate (or Hermitian) transpose.""" 

2646 return self.H 

2647 

2648 @deprecated("2.0", reason="PICOS expressions are now immutable.") 

2649 def copy(self): 

2650 """Return a deep copy of the expression.""" 

2651 from copy import copy as cp 

2652 

2653 return self._basetype(cp(self._symbStr), self._shape, 

2654 {mtbs: cp(coef) for mtbs, coef in self._coefs.items()}) 

2655 

2656 @deprecated("2.0", reason="PICOS expressions are now immutable.") 

2657 def soft_copy(self): 

2658 """Return a shallow copy of the expression.""" 

2659 return self._basetype(self._symbStr, self._shape, self._coefs) 

2660 

2661 

2662# -------------------------------------- 

2663__all__ = api_end(_API_START, globals())