Coverage for picos/expressions/exp_biaffine.py: 82.57%

Shortcuts on this page

r m x   toggle line displays

j k   next/prev highlighted chunk

0   (zero) top of page

1   (one) first highlighted chunk

1205 statements  

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) 

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 as sparse matrices.""" 

367 return { 

368 mtbs: c if isinstance(c, cvxopt.spmatrix) else cvxopt.sparse(c) 

369 for mtbs, c in self._coefs.items()} 

370 

371 @cached_property 

372 def _sparse_linear_coefs(self): 

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

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

375 

376 def _reglyphed_string(self, mutable_name_map): 

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

378 string = self.string 

379 

380 if isinstance(string, glyphs.GlStr): 

381 string = string.reglyphed(mutable_name_map) 

382 elif string in mutable_name_map: 

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

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

385 string = mutable_name_map[string] 

386 

387 return string 

388 

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

390 # Abstract method implementations and overridings for Expression. 

391 # -------------------------------------------------------------------------- 

392 

393 def _get_value(self): 

394 # Create a copy of the constant term. 

395 value = self._constant_coef[:] 

396 

397 # Add value of the linear part. 

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

399 summand = coef * mtb._get_internal_value() 

400 

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

402 value += summand 

403 else: 

404 # Exactly one of the matrices is sparse. 

405 value = value + summand 

406 

407 # Add value of the bilinear part. 

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

409 xValue = x._get_internal_value() 

410 yValue = y._get_internal_value() 

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

412 

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

414 value += summand 

415 else: 

416 # Exactly one of the matrices is sparse. 

417 value = value + summand 

418 

419 # Resize the value to the proper shape. 

420 value.size = self.shape 

421 

422 return value 

423 

424 def _get_shape(self): 

425 return self._shape 

426 

427 @cached_unary_operator 

428 def _get_mutables(self): 

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

430 

431 def _is_convex(self): 

432 return not self._bilinear_coefs 

433 

434 def _is_concave(self): 

435 return not self._bilinear_coefs 

436 

437 def _replace_mutables(self, mapping): 

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

439 if self in mapping: 

440 return mapping[self] 

441 

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

443 string = self._reglyphed_string(name_map) 

444 

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

446 # Fast implementation for the basic case. 

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

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

449 else: 

450 # Turn full mapping into an effective mapping. 

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

452 

453 coefs = {} 

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

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

456 self._update_coefs(coefs, old_mtbs, old_coef) 

457 elif len(old_mtbs) == 1: 

458 assert old_mtbs[0] in mapping 

459 new = mapping[old_mtbs[0]] 

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

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

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

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

464 new = mapping[old_mtbs[0]] 

465 old_mtb = old_mtbs[1] 

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

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

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

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

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

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

472 new = mapping[old_mtbs[1]] 

473 old_mtb = old_mtbs[0] 

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

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

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

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

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

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

480 new1 = mapping[old_mtbs[0]] 

481 new2 = mapping[old_mtbs[1]] 

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

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

484 else: 

485 raise NotImplementedError( 

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

487 "supported unless both are replaced with mutables. " 

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

489 else: 

490 assert False 

491 

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

493 

494 def _freeze_mutables(self, freeze): 

495 string = self._reglyphed_string( 

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

497 

498 coefs = {} 

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

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

501 self._update_coefs(coefs, mtbs, coef) 

502 elif len(mtbs) == 1: 

503 assert mtbs[0] in freeze 

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

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

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

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

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

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

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

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

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

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

514 else: 

515 assert False 

516 

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

518 

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

520 # Python special method implementations. 

521 # -------------------------------------------------------------------------- 

522 

523 def __len__(self): 

524 # Faster version that overrides Expression.__len__. 

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

526 

527 def __getitem__(self, index): 

528 def slice2range(s, length): 

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

530 assert isinstance(s, slice) 

531 

532 # Plug in slice's default values. 

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

534 if ss > 0: 

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

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

537 else: 

538 assert ss < 0 

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

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

541 

542 # Wrap negative indices (once). 

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

544 if sb is None: 

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

546 rb = -1 

547 else: 

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

549 

550 # Clamp out-of-bound indices. 

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

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

553 

554 r = range(ra, rb, ss) 

555 

556 if not r: 

557 raise IndexError("Empty slice.") 

558 

559 return r 

560 

561 def range2slice(r, length): 

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

563 

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

565 """ 

566 assert isinstance(r, range) 

567 

568 if not r: 

569 raise IndexError("Empty range.") 

570 

571 ra = r.start 

572 rb = r.stop 

573 rs = r.step 

574 

575 if rs > 0: 

576 if ra < 0 or rb > length: 

577 raise ValueError( 

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

579 else: 

580 assert rs < 0 

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

582 raise ValueError( 

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

584 

585 if rb == -1: 

586 rb = None 

587 

588 return slice(ra, rb, rs) 

589 

590 def list2slice(l, length): 

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

592 

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

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

595 """ 

596 return range2slice(detect_range(l), length) 

597 

598 def slice2str(s): 

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

600 assert isinstance(s, slice) 

601 

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

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

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

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

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

607 return startStr 

608 else: 

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

610 else: 

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

612 

613 def list2str(l, length): 

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

615 assert isinstance(l, list) 

616 

617 # Extract integers wrapped in a list. 

618 if len(l) == 1: 

619 return str(l[0]) 

620 

621 # Try to convert the list to a slice. 

622 try: 

623 l = list2slice(l, length) 

624 except (ValueError, RuntimeError): 

625 pass 

626 

627 if isinstance(l, list): 

628 if len(l) > 4: 

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

630 else: 

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

632 else: 

633 return slice2str(l) 

634 

635 def any2str(a, length): 

636 if isinstance(a, slice): 

637 return slice2str(a) 

638 elif isinstance(a, list): 

639 return list2str(a, length) 

640 else: 

641 assert False 

642 

643 m, n = self._shape 

644 indexStr = None 

645 isIntList = False 

646 

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

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

649 index = list(index) 

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

651 if len(index) != 2: 

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

653 "exactly two keys.") 

654 

655 try: 

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

657 except TypeError as error: 

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

659 "must be comparable.") from error 

660 

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

662 

663 try: 

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

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

666 

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

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

669 "but a proper matrix.") 

670 

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

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

673 

674 I, J = list(I), list(J) 

675 except Exception as error: 

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

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

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

679 from None 

680 

681 # Represent the selection as a global index list. 

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

683 

684 # Use a special index string. 

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

686 

687 # Don't invoke load_dense_data on the list. 

688 isIntList = True 

689 else: # Global indexing. 

690 index = [index] 

691 

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

693 if not index: 

694 raise IndexError("Empty index.") 

695 elif len(index) > 2: 

696 raise IndexError( 

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

698 

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

700 for axis, selection in enumerate(index): 

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

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

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

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

705 try: 

706 matrix = load_dense_data( 

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

708 

709 if 1 not in matrix.size: 

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

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

712 

713 selection = list(matrix) 

714 except Exception as error: 

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

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

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

718 from None 

719 

720 index[axis] = selection 

721 

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

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

724 index = index[0] 

725 

726 if isinstance(index, slice): 

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

728 else: 

729 shape = len(index) 

730 

731 if indexStr is None: 

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

733 else: # Multiple axis slicing. 

734 if indexStr is None: 

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

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

737 

738 # Convert to a global index list. 

739 RC, shape = [], [] 

740 for axis, selection in enumerate(index): 

741 k = self._shape[axis] 

742 

743 if isinstance(selection, slice): 

744 # Turn the slice into an iterable range. 

745 selection = slice2range(selection, k) 

746 

747 # All indices from a slice are nonnegative. 

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

749 

750 if isinstance(selection, list): 

751 # Wrap once. This is consistent with CVXOPT. 

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

753 

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

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

756 raise IndexError( 

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

758 

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

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

761 raise IndexError( 

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

763 

764 RC.append(selection) 

765 shape.append(len(selection)) 

766 

767 rows, cols = RC 

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

769 

770 # Finalize the string. 

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

772 

773 # Retrieve new coefficients and constant term. 

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

775 

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

777 

778 @convert_operands(sameShape=True) 

779 @refine_operands(stop_at_affine=True) 

780 def __add__(self, other): 

781 if not isinstance(other, BiaffineExpression): 

782 return Expression.__add__(self, other) 

783 

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

785 

786 coefs = {} 

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

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

789 if mtbs in other._coefs else coef 

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

791 coefs.setdefault(mtbs, coef) 

792 

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

794 

795 @convert_operands(sameShape=True) 

796 @refine_operands(stop_at_affine=True) 

797 def __radd__(self, other): 

798 if not isinstance(other, BiaffineExpression): 

799 return Expression.__radd__(self, other) 

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 if not isinstance(other, BiaffineExpression): 

807 return Expression.__sub__(self, other) 

808 

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

810 

811 coefs = {} 

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

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

814 if mtbs in other._coefs else coef 

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

816 coefs.setdefault(mtbs, -coef) 

817 

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

819 

820 @convert_operands(sameShape=True) 

821 @refine_operands(stop_at_affine=True) 

822 def __rsub__(self, other): 

823 if not isinstance(other, BiaffineExpression): 

824 return Expression.__rsub__(self, other) 

825 

826 return other.__sub__(self) 

827 

828 @cached_selfinverse_unary_operator 

829 def __neg__(self): 

830 string = glyphs.clever_neg(self.string) 

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

832 

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

834 

835 @convert_operands(rMatMul=True) 

836 @refine_operands(stop_at_affine=True) 

837 def __mul__(self, other): 

838 if not isinstance(other, BiaffineExpression): 

839 return Expression.__mul__(self, other) 

840 

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

842 

843 m, n = self._shape 

844 p, q = other._shape 

845 

846 # Handle or prepare multiplication with a scalar. 

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

848 if self.constant or other.constant: 

849 # Handle all cases involving a constant operand immediately. 

850 if self.constant: 

851 lhs = self._constant_coef 

852 RHS = other._coefs 

853 else: 

854 lhs = other._constant_coef 

855 RHS = self._coefs 

856 

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

858 

859 # Work around CVXOPT not broadcasting sparse scalars. 

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

861 lhs = lhs[0] 

862 

863 return self._common_basetype(other)( 

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

865 elif n == p: 

866 # Regular matrix multiplication already works. 

867 pass 

868 elif m*n == 1: 

869 # Prepare a regular matrix multiplication. 

870 # HACK: Replaces self. 

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

872 m, n = self._shape 

873 else: 

874 # Prepare a regular matrix multiplication. 

875 assert p*q == 1 

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

877 p, q = other._shape 

878 

879 assert n == p 

880 

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

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

883 row_by_column = m == 1 and q == 1 

884 

885 # Regular matrix multiplication. 

886 coefs = {} 

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

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

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

890 raise NotImplementedError("Product not supported: " 

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

892 elif len(lmtbs) == 0: 

893 # Constant by constant, linear or bilinear. 

894 if row_by_column: 

895 factor = lcoef.T 

896 else: 

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

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

899 

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

901 elif len(rmtbs) == 0: 

902 # Constant, linear or bilinear by constant. 

903 if row_by_column: 

904 factor = rcoef.T 

905 else: 

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

907 factor = right_kronecker_I( 

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

909 

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

911 else: 

912 # Linear by linear. 

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

914 

915 if row_by_column: 

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

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

918 else: 

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

920 # (Stahlberg 2020). 

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

922 

923 A = cvxopt_K(m, n)*lcoef 

924 B = rcoef 

925 

926 A = A.T 

927 B = B.T 

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

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

930 

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

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

933 A.size = n, m*a 

934 B.size = p, q*b 

935 A = A.T 

936 

937 C = A*B 

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

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

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

941 C = C.T 

942 

943 coef = C 

944 

945 # Forbid quadratic results. 

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

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

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

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

950 

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

952 

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

954 

955 @convert_operands(lMatMul=True) 

956 @refine_operands(stop_at_affine=True) 

957 def __rmul__(self, other): 

958 if not isinstance(other, BiaffineExpression): 

959 return Expression.__rmul__(self, other) 

960 

961 return other.__mul__(self) 

962 

963 @convert_operands(sameShape=True) 

964 @refine_operands(stop_at_affine=True) 

965 def __or__(self, other): 

966 if not isinstance(other, BiaffineExpression): 

967 return Expression.__or__(self, other) 

968 

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

970 result._symbStr = glyphs.clever_dotp( 

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

972 return result 

973 

974 @convert_operands(sameShape=True) 

975 @refine_operands(stop_at_affine=True) 

976 def __ror__(self, other): 

977 if not isinstance(other, BiaffineExpression): 

978 return Expression.__ror__(self, other) 

979 

980 return other.__or__(self) 

981 

982 @convert_operands(sameShape=True) 

983 @refine_operands(stop_at_affine=True) 

984 def __xor__(self, other): 

985 if not isinstance(other, BiaffineExpression): 

986 return Expression.__xor__(self, other) 

987 

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

989 

990 if other.constant: 

991 factor = cvxopt.spdiag(other._constant_coef) 

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

993 elif self.constant: 

994 factor = cvxopt.spdiag(self._constant_coef) 

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

996 else: 

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

998 

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

1000 

1001 @convert_operands(sameShape=True) 

1002 @refine_operands(stop_at_affine=True) 

1003 def __rxor__(self, other): 

1004 if not isinstance(other, BiaffineExpression): 

1005 return Expression.__rxor__(self, other) 

1006 

1007 return other.__xor__(self) 

1008 

1009 @convert_operands() 

1010 @refine_operands(stop_at_affine=True) 

1011 def __matmul__(self, other): 

1012 if not isinstance(other, BiaffineExpression): 

1013 return Expression.__matmul__(self, other) 

1014 

1015 cls = self._common_basetype(other) 

1016 tc = cls._get_typecode() 

1017 

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

1019 

1020 m, n = self._shape 

1021 p, q = other._shape 

1022 

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

1024 Kqm = cvxopt_K(q, m) 

1025 Kqm_Ip = right_kronecker_I(Kqm, p) 

1026 In_Kqm_Ip = left_kronecker_I(Kqm_Ip, n) 

1027 

1028 def _kron(A, B): 

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

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

1031 return In_Kqm_Ip*A_B*Kab 

1032 

1033 coefs = {} 

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

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

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

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

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

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

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

1041 

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

1043 else: 

1044 raise NotImplementedError("Product not supported: " 

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

1046 

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

1048 

1049 @convert_operands() 

1050 @refine_operands(stop_at_affine=True) 

1051 def __rmatmul__(self, other): 

1052 if not isinstance(other, BiaffineExpression): 

1053 return Expression.__rmatmul__(self, other) 

1054 

1055 return other.__matmul__(self) 

1056 

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

1058 def kron(self, other): 

1059 """Denote the Kronecker product with another expression on the right.""" 

1060 return self.__matmul__(other) 

1061 

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

1063 def leftkron(self, other): 

1064 """Denote the Kronecker product with another expression on the left.""" 

1065 return self.__rmatmul__(other) 

1066 

1067 @convert_operands(scalarRHS=True) 

1068 @refine_operands(stop_at_affine=True) 

1069 def __truediv__(self, other): 

1070 if not isinstance(other, BiaffineExpression): 

1071 return Expression.__truediv__(self, other) 

1072 

1073 if not other.constant: 

1074 raise TypeError( 

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

1076 

1077 if other.is0: 

1078 raise ZeroDivisionError( 

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

1080 

1081 divisor = other._constant_coef[0] 

1082 

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

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

1085 

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

1087 

1088 @convert_operands(scalarLHS=True) 

1089 @refine_operands(stop_at_affine=True) 

1090 def __rtruediv__(self, other): 

1091 if not isinstance(other, BiaffineExpression): 

1092 return Expression.__rtruediv__(self, other) 

1093 

1094 return other.__truediv__(self) 

1095 

1096 @convert_operands(horiCat=True) 

1097 @refine_operands(stop_at_affine=True) 

1098 def __and__(self, other): 

1099 if not isinstance(other, BiaffineExpression): 

1100 return Expression.__and__(self, other) 

1101 

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

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

1104 

1105 coefs = {} 

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

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

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

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

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

1111 

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

1113 

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

1115 

1116 @convert_operands(horiCat=True) 

1117 @refine_operands(stop_at_affine=True) 

1118 def __rand__(self, other): 

1119 if not isinstance(other, BiaffineExpression): 

1120 return Expression.__rand__(self, other) 

1121 

1122 return other.__and__(self) 

1123 

1124 @convert_operands(vertCat=True) 

1125 @refine_operands(stop_at_affine=True) 

1126 def __floordiv__(self, other): 

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

1128 p, q = upperRows, lowerRows 

1129 return [column for columnPairs in [ 

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

1131 for column in columnPairs] 

1132 

1133 if not isinstance(other, BiaffineExpression): 

1134 return Expression.__floordiv__(self, other) 

1135 

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

1137 

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

1139 shape = (p + q, c) 

1140 

1141 coefs = {} 

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

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

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

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

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

1147 

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

1149 

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

1151 

1152 @convert_operands(vertCat=True) 

1153 @refine_operands(stop_at_affine=True) 

1154 def __rfloordiv__(self, other): 

1155 if not isinstance(other, BiaffineExpression): 

1156 return Expression.__rfloordiv__(self, other) 

1157 

1158 return other.__floordiv__(self) 

1159 

1160 @convert_operands(scalarRHS=True) 

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

1162 def __pow__(self, other): 

1163 from .exp_powtrace import PowerTrace 

1164 

1165 if not self.scalar: 

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

1167 

1168 if not other.constant: 

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

1170 

1171 if other.complex: 

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

1173 

1174 exponent = other.value 

1175 

1176 if exponent == 2: 

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

1178 else: 

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

1180 

1181 # -------------------------------------------------------------------------- 

1182 # Properties and functions that describe the expression. 

1183 # -------------------------------------------------------------------------- 

1184 

1185 @cached_property 

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

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

1188 

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

1190 

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

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

1193 larger numeric errors or the effects of noisy data. 

1194 """ 

1195 return self.equals( 

1196 self.H, relTol=settings.RELATIVE_HERMITIANNESS_TOLERANCE) 

1197 

1198 @property 

1199 def is0(self): 

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

1201 return not self._coefs 

1202 

1203 @cached_property 

1204 def is1(self): 

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

1206 # Must be a scalar or vector. 

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

1208 return False 

1209 

1210 # Must be constant. 

1211 if not self.constant: 

1212 return False 

1213 

1214 # Constant term must be all ones. 

1215 return not self._constant_coef - 1 

1216 

1217 @cached_property 

1218 def isI(self): 

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

1220 m, n = self._shape 

1221 

1222 # Must be a square matrix. 

1223 if m != n: 

1224 return False 

1225 

1226 # Must be constant. 

1227 if not self.constant: 

1228 return False 

1229 

1230 # Constant term must be the identity. 

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

1232 

1233 @cached_property 

1234 def complex(self): 

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

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

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

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

1239 

1240 @property 

1241 def isreal(self): 

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

1243 return not self.complex 

1244 

1245 @convert_operands() 

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

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

1248 

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

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

1251 coefficients and constant term can be equal to an 

1252 :class:`~.exp_affine.AffineExpression`. 

1253 

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

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

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

1257 this expression. In particular, 

1258 

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

1260 - lists and tuples are treated as column vectors and 

1261 - algebraic strings must specify a shape (see 

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

1263 

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

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

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

1267 entries of the coefficient matrices and the constant terms being 

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

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

1270 entries of the coefficient matrices and the constant terms being 

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

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

1273 

1274 :Example: 

1275 

1276 >>> from picos import Constant 

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

1278 >>> repr(A) 

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

1280 >>> A.is0 

1281 True 

1282 >>> A.equals(0) 

1283 False 

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

1285 True 

1286 >>> repr(A*1j) 

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

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

1289 True 

1290 """ 

1291 if self is other: 

1292 return True 

1293 

1294 if not isinstance(other, BiaffineExpression): 

1295 return False 

1296 

1297 if self._shape != other._shape: 

1298 return False 

1299 

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

1301 for mtbs in self._bilinear_coefs), \ 

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

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

1304 

1305 for mtbs in other._coefs: 

1306 if mtbs not in self._coefs: 

1307 return False 

1308 

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

1310 if mtbs not in other._coefs: 

1311 return False 

1312 

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

1314 return False 

1315 

1316 return True 

1317 

1318 # -------------------------------------------------------------------------- 

1319 # Methods and properties that return modified copies. 

1320 # -------------------------------------------------------------------------- 

1321 

1322 def renamed(self, string): 

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

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

1325 

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

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

1328 

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

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

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

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

1333 generalization of row-major) by default. 

1334 

1335 :param str order: 

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

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

1338 

1339 :Example: 

1340 

1341 >>> from picos import Constant 

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

1343 >>> print(C) 

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

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

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

1347 [ 0.00e+00 3.00e+00] 

1348 [ 1.00e+00 4.00e+00] 

1349 [ 2.00e+00 5.00e+00] 

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

1351 [ 0.00e+00 2.00e+00] 

1352 [ 4.00e+00 1.00e+00] 

1353 [ 3.00e+00 5.00e+00] 

1354 """ 

1355 if order not in "FC": 

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

1357 

1358 shape = load_shape(shape, wildcards=True) 

1359 

1360 if shape == self._shape: 

1361 return self 

1362 elif shape == (None, None): 

1363 return self 

1364 

1365 length = len(self) 

1366 

1367 if shape[0] is None: 

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

1369 elif shape[1] is None: 

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

1371 

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

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

1374 

1375 if order == "F": 

1376 coefs = self._coefs 

1377 reshaped_glyph = glyphs.reshaped 

1378 else: 

1379 m, n = self._shape 

1380 p, q = shape 

1381 

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

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

1384 R = K_new * K_old 

1385 

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

1387 reshaped_glyph = glyphs.reshaprm 

1388 

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

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

1391 

1392 def broadcasted(self, shape): 

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

1394 

1395 :Example: 

1396 

1397 >>> from picos import Constant 

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

1399 >>> print(C) 

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

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

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

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

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

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

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

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

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

1409 """ 

1410 shape = load_shape(shape, wildcards=True) 

1411 shape = blend_shapes(shape, self._shape) 

1412 

1413 if shape == self._shape: 

1414 return self 

1415 

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

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

1418 

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

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

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

1422 

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

1424 string = glyphs.matrix(self.string) 

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

1426 

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

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

1429 

1430 def reshaped_or_broadcasted(self, shape): 

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

1432 

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

1434 may not contain a wildcard character. 

1435 

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

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

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

1439 """ 

1440 shape = load_shape(shape, wildcards=False) 

1441 

1442 try: 

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

1444 return self.reshaped(shape) 

1445 else: 

1446 return self.broadcasted(shape) 

1447 except ValueError: 

1448 raise ValueError( 

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

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

1451 

1452 @cached_property 

1453 def hermitianized(self): 

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

1455 

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

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

1458 

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

1460 subspace of symmetric matrices. 

1461 """ 

1462 if not self.square: 

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

1464 

1465 return (self + self.H)/2 

1466 

1467 @cached_property 

1468 def real(self): 

1469 """Real part of the expression.""" 

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

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

1472 

1473 @cached_property 

1474 def imag(self): 

1475 """Imaginary part of the expression.""" 

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

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

1478 

1479 @cached_property 

1480 def bilin(self): 

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

1482 return self._basetype( 

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

1484 

1485 @cached_property 

1486 def lin(self): 

1487 """Linear part of the expression.""" 

1488 return self._basetype( 

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

1490 

1491 @cached_property 

1492 def noncst(self): 

1493 """Nonconstant part of the expression.""" 

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

1495 

1496 return self._basetype( 

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

1498 

1499 @cached_property 

1500 def cst(self): 

1501 """Constant part of the expression.""" 

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

1503 

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

1505 

1506 @cached_selfinverse_property 

1507 def T(self): 

1508 """Matrix transpose.""" 

1509 if len(self) == 1: 

1510 return self 

1511 

1512 m, n = self._shape 

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

1514 

1515 string = glyphs.transp(self.string) 

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

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

1518 

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

1520 

1521 @cached_selfinverse_property 

1522 def conj(self): 

1523 """Complex conjugate.""" 

1524 string = glyphs.conj(self.string) 

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

1526 

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

1528 

1529 @cached_selfinverse_property 

1530 def H(self): 

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

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

1533 

1534 def _square_equal_subsystem_dims(self, diagLen): 

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

1536 m, n = self._shape 

1537 k = math.log(m, diagLen) 

1538 

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

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

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

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

1543 

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

1545 

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

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

1548 

1549 If the expression can be written as 

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

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

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

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

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

1555 

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

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

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

1559 :type subsystems: int or tuple or list 

1560 

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

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

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

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

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

1566 :type dimensions: int or tuple or list 

1567 

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

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

1570 

1571 :Example: 

1572 

1573 >>> from picos import Constant 

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

1594 """ 

1595 m, n = self._shape 

1596 

1597 if isinstance(dimensions, int): 

1598 dimensions = self._square_equal_subsystem_dims(dimensions) 

1599 else: 

1600 dimensions = [ 

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

1602 

1603 if reduce( 

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

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

1606 

1607 if isinstance(subsystems, int): 

1608 subsystems = (subsystems,) 

1609 elif not subsystems: 

1610 return self 

1611 

1612 numSys = len(dimensions) 

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

1614 

1615 for sys in subsystems: 

1616 if not isinstance(sys, int): 

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

1618 .format(type(sys).__name__)) 

1619 elif sys < 0: 

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

1621 elif sys >= numSys: 

1622 raise IndexError( 

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

1624 .format(sys, numSys)) 

1625 

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

1627 if len(subsystems) == numSys: 

1628 return self.T 

1629 

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

1631 d = m * n 

1632 V = [1]*d 

1633 I = range(d) 

1634 J = cvxopt.matrix(I) 

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

1636 obh, obw = 1, 1 

1637 sysStrings = None 

1638 

1639 # Apply transpositions starting with the rightmost Kronecker factor. 

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

1641 # Shape of current system. 

1642 p, q = dimensions[sys] 

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

1644 

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

1646 ibh, ibw = obh, obw 

1647 

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

1649 # maintained but that are subject to transposition independently. 

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

1651 obh *= p 

1652 obw *= q 

1653 

1654 # Only transpose selected subsystems. 

1655 if sys not in subsystems: 

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

1657 if sysStrings else sysString 

1658 continue 

1659 else: 

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

1661 if sysStrings else glyphs.transp(sysString) 

1662 

1663 # Shape of outer blocks after transposition. 

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

1665 

1666 # Shape of full matrix after transposition. 

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

1668 

1669 for vi in I: 

1670 # Full matrix column and row. 

1671 c, r = divmod(vi, m) 

1672 

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

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

1675 obi, obr = divmod(r, obh) 

1676 obj, obc = divmod(c, obw) 

1677 

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

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

1680 ibi, ibr = divmod(obr, ibh) 

1681 ibj, ibc = divmod(obc, ibw) 

1682 

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

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

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

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

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

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

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

1690 

1691 # Prepare the transposition. 

1692 T[tvi] = J[vi] 

1693 

1694 # Apply the transposition. 

1695 J, T = T, J 

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

1697 

1698 # Finalize the partial commutation matrix. 

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

1700 

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

1702 shape = (m, n) 

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

1704 

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

1706 

1707 @cached_property 

1708 def T0(self): 

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

1710 

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

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

1713 """ 

1714 return self.partial_transpose(subsystems=0) 

1715 

1716 @cached_property 

1717 def T1(self): 

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

1719 

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

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

1722 """ 

1723 return self.partial_transpose(subsystems=1) 

1724 

1725 @cached_property 

1726 def T2(self): 

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

1728 

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

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

1731 """ 

1732 return self.partial_transpose(subsystems=2) 

1733 

1734 @cached_property 

1735 def T3(self): 

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

1737 

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

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

1740 """ 

1741 return self.partial_transpose(subsystems=3) 

1742 

1743 @cached_property 

1744 def Tl(self): 

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

1746 

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

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

1749 """ 

1750 return self.partial_transpose(subsystems=-1) 

1751 

1752 @staticmethod 

1753 def _reindex_F(indices, source, destination): 

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

1755 new = [] 

1756 offset = 0 

1757 factor = 1 

1758 

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

1760 offset += factor*index 

1761 factor *= base 

1762 

1763 for base in destination: 

1764 offset, remainder = divmod(offset, base) 

1765 new.append(remainder) 

1766 

1767 return tuple(new) 

1768 

1769 @staticmethod 

1770 def _reindex_C(indices, source, destination): 

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

1772 new = [] 

1773 offset = 0 

1774 factor = 1 

1775 

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

1777 offset += factor*index 

1778 factor *= base 

1779 

1780 for base in reversed(destination): 

1781 offset, remainder = divmod(offset, base) 

1782 new.insert(0, remainder) 

1783 

1784 return tuple(new) 

1785 

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

1787 """Return the reshuffled or realigned expression. 

1788 

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

1790 the following sequence of operations: 

1791 

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

1793 according to ``order``. 

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

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

1796 according to ``order``. 

1797 

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

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

1800 

1801 .. code:: 

1802 

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

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

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

1806 reshuffled = numpy.einsum(P, reshuffled) 

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

1808 

1809 :param permutation: 

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

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

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

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

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

1815 should be swapped. 

1816 :type permutation: 

1817 str or tuple or list 

1818 

1819 :param dimensions: 

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

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

1822 hypercubic and the number of dimensions is inferred from the 

1823 ``permutation`` argument. 

1824 :type dimensions: 

1825 None or tuple or list 

1826 

1827 :param str order: 

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

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

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

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

1832 

1833 :Example: 

1834 

1835 >>> from picos import Constant 

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

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

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

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

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

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

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

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

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

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

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

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

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

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

1850 True 

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

1852 True 

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

1854 True 

1855 """ 

1856 m, n = self._shape 

1857 mn = m*n 

1858 

1859 # Load the permutation. 

1860 ordered = sorted(permutation) 

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

1862 

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

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

1865 "contain duplicate elements.") 

1866 

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

1868 

1869 numDims = len(P) 

1870 

1871 # Load the dimensions. 

1872 guessDimensions = dimensions is None 

1873 

1874 if guessDimensions: 

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

1876 else: 

1877 if len(dimensions) != numDims: 

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

1879 "number of dimensions.") 

1880 

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

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

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

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

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

1886 

1887 # Load the indexing order. 

1888 if order not in "FC": 

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

1890 

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

1892 

1893 # Nothing to do for the neutral permutation. 

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

1895 return self 

1896 

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

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

1899 for i in range(mn): 

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

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

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

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

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

1905 I.append(k*m + j) 

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

1907 

1908 # Create the string. 

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

1910 

1911 if not guessDimensions: 

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

1913 

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

1915 

1916 # Finalize the new expression. 

1917 shape = (m, n) 

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

1919 

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

1921 

1922 @cached_property 

1923 def sum(self): 

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

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

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

1927 return (self | 1) 

1928 

1929 @cached_property 

1930 def rowsum(self): 

1931 """Sum over the rows of the expression as a row vector.""" 

1932 from .algebra import J 

1933 

1934 return J(self.shape[0]).T * self 

1935 

1936 @cached_property 

1937 def colsum(self): 

1938 """Sum over the columns of the expression as a column vector.""" 

1939 from .algebra import J 

1940 

1941 return self * J(self.shape[1]) 

1942 

1943 @cached_property 

1944 def tr(self): 

1945 """Trace of a square expression.""" 

1946 if not self.square: 

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

1948 .format(self.string)) 

1949 

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

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

1952 return (self | "I") 

1953 

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

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

1956 

1957 If the expression can be written as 

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

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

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

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

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

1963 

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

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

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

1967 :type subsystems: int or tuple or list 

1968 

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

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

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

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

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

1974 :type dimensions: int or tuple or list 

1975 

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

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

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

1979 way. 

1980 

1981 :Example: 

1982 

1983 >>> from picos import Constant 

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

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

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

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

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

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

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

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

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

1993 [ 1.00e+01 1.80e+01] 

1994 [ 1.20e+01 2.00e+01] 

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

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

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

1998 [ 5.00e+00 2.10e+01] 

1999 [ 9.00e+00 2.50e+01] 

2000 """ 

2001 # Shape of the original matrix. 

2002 m, n = self._shape 

2003 

2004 if isinstance(dimensions, int): 

2005 dimensions = self._square_equal_subsystem_dims(dimensions) 

2006 else: 

2007 dimensions = [ 

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

2009 

2010 if reduce( 

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

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

2013 

2014 if isinstance(subsystems, int): 

2015 subsystems = (subsystems,) 

2016 elif not subsystems: 

2017 return self 

2018 

2019 numSys = len(dimensions) 

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

2021 

2022 for sys in subsystems: 

2023 if not isinstance(sys, int): 

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

2025 .format(type(sys).__name__)) 

2026 elif sys < 0: 

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

2028 elif sys >= numSys: 

2029 raise IndexError( 

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

2031 .format(sys, numSys)) 

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

2033 raise TypeError( 

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

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

2036 

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

2038 if len(subsystems) == numSys: 

2039 return self.tr 

2040 

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

2042 T = [] 

2043 

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

2045 ibh, ibw = m, n 

2046 sysStrings = None 

2047 for sys in range(numSys): 

2048 # Shape of current system. 

2049 p, q = dimensions[sys] 

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

2051 

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

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

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

2055 obh, obw = ibh, ibw 

2056 

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

2058 # main-diagonal blocks of an outer block. 

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

2060 

2061 # Only trace over selected subsystems. 

2062 if sys not in subsystems: 

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

2064 if sysStrings else sysString 

2065 continue 

2066 else: 

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

2068 if sysStrings else glyphs.trace(sysString) 

2069 

2070 # Shape of new matrix. 

2071 assert p == q 

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

2073 

2074 # Prepare one factor of T. 

2075 oldLen = m * n 

2076 newLen = mN * nN 

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

2078 shape = (newLen, oldLen) 

2079 

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

2081 for viN in range(newLen): 

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

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

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

2085 cN, rN = divmod(viN, mN) 

2086 

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

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

2089 obi, ibr = divmod(rN, ibh) 

2090 obj, ibc = divmod(cN, ibw) 

2091 

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

2093 # scalar per on-diagonal inner block. 

2094 for k in range(p): 

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

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

2097 I.append(viN) 

2098 J.append(cO*m + rO) 

2099 

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

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

2102 

2103 # Make next iteration work on the new matrix. 

2104 m, n = mN, nN 

2105 

2106 # Multiply out the linear partial trace operator T. 

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

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

2109 

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

2111 shape = (m, n) 

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

2113 

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

2115 

2116 @cached_property 

2117 def tr0(self): 

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

2119 

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

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

2122 """ 

2123 return self.partial_trace(subsystems=0) 

2124 

2125 @cached_property 

2126 def tr1(self): 

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

2128 

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

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

2131 """ 

2132 return self.partial_trace(subsystems=1) 

2133 

2134 @cached_property 

2135 def tr2(self): 

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

2137 

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

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

2140 """ 

2141 return self.partial_trace(subsystems=2) 

2142 

2143 @cached_property 

2144 def tr3(self): 

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

2146 

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

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

2149 """ 

2150 return self.partial_trace(subsystems=3) 

2151 

2152 @cached_property 

2153 def trl(self): 

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

2155 

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

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

2158 """ 

2159 return self.partial_trace(subsystems=-1) 

2160 

2161 @cached_property 

2162 def vec(self): 

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

2164 

2165 .. note:: 

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

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

2168 its result is cached. 

2169 

2170 :Example: 

2171 

2172 >>> from picos import Constant 

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

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

2175 True 

2176 >>> A[:] is A[:] 

2177 False 

2178 >>> A.vec is A.vec 

2179 True 

2180 """ 

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

2182 return self 

2183 else: 

2184 return self._basetype( 

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

2186 

2187 def dupvec(self, n): 

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

2189 

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

2191 

2192 :returns: A column vector. 

2193 

2194 :Example: 

2195 

2196 >>> from picos import Constant 

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

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

2199 True 

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

2201 True 

2202 """ 

2203 if not isinstance(n, int): 

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

2205 

2206 if n < 1: 

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

2208 

2209 if n == 1: 

2210 return self.vec 

2211 else: 

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

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

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

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

2216 

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

2218 

2219 @cached_property 

2220 def trilvec(self): 

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

2222 

2223 :returns: 

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

2225 :math:`i \geq j`. 

2226 

2227 .. note:: 

2228 

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

2230 instead of ``A.trilvec``. 

2231 

2232 :Example: 

2233 

2234 >>> from picos import Constant 

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

2236 >>> print(A) 

2237 [ 1.00e+00 2.00e+00] 

2238 [ 3.00e+00 4.00e+00] 

2239 [ 5.00e+00 6.00e+00] 

2240 >>> print(A.trilvec) 

2241 [ 1.00e+00] 

2242 [ 3.00e+00] 

2243 [ 5.00e+00] 

2244 [ 4.00e+00] 

2245 [ 6.00e+00] 

2246 """ 

2247 m, n = self._shape 

2248 

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

2250 return self 

2251 elif m == 1: # Row vector. 

2252 return self[0] 

2253 

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

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

2256 d = len(rows) 

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

2258 

2259 string = glyphs.trilvec(self.string) 

2260 shape = (d, 1) 

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

2262 

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

2264 

2265 @cached_property 

2266 def triuvec(self): 

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

2268 

2269 :returns: 

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

2271 :math:`i \leq j`. 

2272 

2273 .. note:: 

2274 

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

2276 instead of ``A.triuvec``. 

2277 

2278 :Example: 

2279 

2280 >>> from picos import Constant 

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

2282 >>> print(A) 

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

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

2285 >>> print(A.triuvec) 

2286 [ 1.00e+00] 

2287 [ 2.00e+00] 

2288 [ 5.00e+00] 

2289 [ 3.00e+00] 

2290 [ 6.00e+00] 

2291 """ 

2292 m, n = self._shape 

2293 

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

2295 return self 

2296 elif n == 1: # Column vector. 

2297 return self[0] 

2298 

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

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

2301 d = len(rows) 

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

2303 

2304 string = glyphs.triuvec(self.string) 

2305 shape = (d, 1) 

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

2307 

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

2309 

2310 @cached_property 

2311 def svec(self): 

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

2313 

2314 In the real symmetric case 

2315 

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

2317 - the vectorization is isometric and isomorphic, and 

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

2319 :class:`~.variables.SymmetricVariable` class. 

2320 

2321 In the complex hermitian case 

2322 

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

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

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

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

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

2328 

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

2330 

2331 :raises TypeError: 

2332 If the expression is not square. 

2333 

2334 :raises ValueError: 

2335 If the expression is not hermitian. 

2336 """ 

2337 if not self.square: 

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

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

2340 elif not self.hermitian: 

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

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

2343 

2344 vec = SymmetricVectorization(self._shape) 

2345 V = vec._full2special 

2346 

2347 string = glyphs.svec(self.string) 

2348 

2349 if self.isreal: 

2350 vec = SymmetricVectorization(self._shape) 

2351 V = vec._full2special 

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

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

2354 else: 

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

2356 # reasons, SymmetricVectorization averages the elements in the 

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

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

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

2360 real_svec = self.real.svec 

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

2362 

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

2364 

2365 with unlocked_cached_properties(result): 

2366 result.desvec = self 

2367 

2368 return result 

2369 

2370 @cached_property 

2371 def desvec(self): 

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

2373 

2374 :raises TypeError: 

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

2376 permissible set 

2377 

2378 .. math:: 

2379 

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

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

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

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

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

2385 

2386 :raises ValueError: 

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

2388 of :attr:`svec`. 

2389 """ 

2390 if 1 not in self.shape: 

2391 raise TypeError( 

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

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

2394 

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

2396 if int(n) != n: 

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

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

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

2400 n = int(n) 

2401 

2402 vec = SymmetricVectorization((n, n)) 

2403 D = vec._special2full 

2404 string = glyphs.desvec(self.string) 

2405 

2406 if self.isreal: 

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

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

2409 

2410 assert result.hermitian 

2411 else: 

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

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

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

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

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

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

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

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

2420 

2421 real_desvec = self.real.desvec 

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

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

2424 

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

2426 

2427 if not result.hermitian: 

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

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

2430 .format(self.string)) 

2431 

2432 with unlocked_cached_properties(result): 

2433 result.svec = self 

2434 

2435 return result 

2436 

2437 def dupdiag(self, n): 

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

2439 

2440 Vectorization is performed in column-major order. 

2441 

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

2443 """ 

2444 from .algebra import I 

2445 

2446 if self.scalar: 

2447 return self * I(n) 

2448 

2449 # Vectorize and duplicate the expression. 

2450 vec = self.dupvec(n) 

2451 d = len(vec) 

2452 

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

2454 ones = [1]*d 

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

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

2457 

2458 string = glyphs.diag(vec.string) 

2459 shape = (d, d) 

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

2461 

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

2463 

2464 @cached_property 

2465 def diag(self): 

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

2467 

2468 Vectorization is performed in column-major order. 

2469 """ 

2470 from .algebra import O, I 

2471 

2472 if self.is0: 

2473 return O(len(self), len(self)) 

2474 elif self.is1: 

2475 return I(len(self)) 

2476 else: 

2477 return self.dupdiag(1) 

2478 

2479 @cached_property 

2480 def maindiag(self): 

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

2482 if 1 in self._shape: 

2483 return self[0] 

2484 

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

2486 step = self._shape[0] + 1 

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

2488 d = len(rows) 

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

2490 

2491 string = glyphs.maindiag(self.string) 

2492 shape = (d, 1) 

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

2494 

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

2496 

2497 def factor_out(self, mutable): 

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

2499 

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

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

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

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

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

2505 

2506 :returns: 

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

2508 on ``mutable``. 

2509 

2510 :raises TypeError: 

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

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

2513 format. 

2514 

2515 :raises LookupError: 

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

2517 

2518 :Example: 

2519 

2520 >>> from picos import RealVariable 

2521 >>> from picos.uncertain import UnitBallPerturbationSet 

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

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

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

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

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

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

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

2529 >>> az 

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

2531 >>> a0 

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

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

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

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

2536 True 

2537 """ 

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

2539 raise TypeError( 

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

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

2542 

2543 mtb = mutable 

2544 

2545 if not isinstance(mtb, Mutable): 

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

2547 .format(type(mtb).__name__)) 

2548 

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

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

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

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

2553 

2554 if mtb not in self.mutables: 

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

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

2557 

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

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

2560 

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

2562 ax_coefs = {} 

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

2564 if mtb not in set(mtbs): 

2565 continue 

2566 elif len(mtbs) == 1: 

2567 assert mtbs[0] is mtb 

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

2569 else: 

2570 if mtbs[0] is mtb: 

2571 other_mtb = mtbs[1] 

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

2573 else: 

2574 assert mtbs[1] is mtb 

2575 other_mtb = mtbs[0] 

2576 C = coef 

2577 

2578 d = other_mtb.dim 

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

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

2581 

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

2583 

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

2585 a0_shape = self._shape 

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

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

2588 

2589 return ax.refined, a0.refined 

2590 

2591 # -------------------------------------------------------------------------- 

2592 # Backwards compatibility methods. 

2593 # -------------------------------------------------------------------------- 

2594 

2595 @classmethod 

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

2597 def fromScalar(cls, scalar): 

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

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

2600 

2601 @classmethod 

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

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

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

2605 return cls.from_constant(matrix, size) 

2606 

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

2608 def hadamard(self, fact): 

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

2610 return self.__xor__(fact) 

2611 

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

2613 def isconstant(self): 

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

2615 return self.constant 

2616 

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

2618 def same_as(self, other): 

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

2620 return self.equals(other) 

2621 

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

2623 def transpose(self): 

2624 """Return the matrix transpose.""" 

2625 return self.T 

2626 

2627 @cached_property 

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

2629 def Tx(self): 

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

2631 m, n = self._shape 

2632 dims = None 

2633 

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

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

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

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

2638 break 

2639 

2640 if dims: 

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

2642 else: 

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

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

2645 .format(glyphs.shape( 

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

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

2648 

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

2650 def conjugate(self): 

2651 """Return the complex conjugate.""" 

2652 return self.conj 

2653 

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

2655 def Htranspose(self): 

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

2657 return self.H 

2658 

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

2660 def copy(self): 

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

2662 from copy import copy as cp 

2663 

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

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

2666 

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

2668 def soft_copy(self): 

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

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

2671 

2672 

2673# -------------------------------------- 

2674__all__ = api_end(_API_START, globals())