Coverage for picos/expressions/exp_biaffine.py: 82.94%
1219 statements
« prev ^ index » next coverage.py v7.6.12, created at 2025-04-12 07:53 +0000
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# ------------------------------------------------------------------------------
20"""Implements :class:`BiaffineExpression`."""
22import math
23import operator
24from abc import ABC, abstractmethod
25from functools import reduce
26from itertools import product
28import cvxopt
29import numpy
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
45_API_START = api_start(globals())
46# -------------------------------
49class BiaffineExpression(Expression, ABC):
50 r"""A multidimensional (complex) biaffine expression.
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.
58 In general this expression has the form
60 .. math::
62 A(x,y) = B(x,y) + P(x) + Q(y) + C
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.
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`.
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 """
86 # --------------------------------------------------------------------------
87 # Static methods.
88 # --------------------------------------------------------------------------
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
99 # --------------------------------------------------------------------------
100 # Initialization and factory methods.
101 # --------------------------------------------------------------------------
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 {}"
108 def __init__(self, string, shape=(1, 1), coefficients={}):
109 """Initialize a (complex) biaffine expression.
111 :param str string: A symbolic string description.
113 :param shape: Shape of a vector or matrix expression.
114 :type shape: int or tuple or list
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.
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).
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
135 shape = load_shape(shape)
136 length = shape[0]*shape[1]
138 if not isinstance(coefficients, dict):
139 raise TypeError("Coefficients of {} must be given as a dict."
140 .format(type(self).__name__))
142 # Store shape.
143 self._shape = shape
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__))
154 if not all(isinstance(mtb, Mutable) for mtb in mtbs):
155 raise TypeError("Coefficients must be indexed by mutables.")
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__))
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__))
174 # Do not store coefficients of zero.
175 if not coef:
176 continue
178 # Dimension of the tensor product of the vectorized mutables.
179 dim = reduce(operator.mul, (mtb.dim for mtb in mtbs), 1)
181 # Load the coefficient matrix in the desired format.
182 coef = load_data(
183 coef, (length, dim), self._typecode, alwaysCopy=False)[0]
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)
190 # Make coef apply to vec(y*x.T) instead of vec(x*y.T).
191 coef = coef * K
193 # Swap x and y.
194 mtbs = (mtbs[1], mtbs[0])
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
202 # Determine the type string.
203 typeStr = self._get_type_string_base()
205 if "{}" in typeStr:
206 hasBilinear = bool(self._bilinear_coefs)
207 hasLinear = bool(self._linear_coefs)
208 hasConstant = bool(self._constant_coef)
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")
223 typeStr = "{} {}".format(glyphs.shape(shape), typeStr)
225 Expression.__init__(self, typeStr, string)
227 @classmethod
228 def from_constant(cls, constant, shape=None, name=None):
229 """Create a class instance from the given numeric constant.
231 Loads the given constant as a PICOS expression, optionally broadcasted
232 or reshaped to the given shape and named as specified.
234 See :func:`~.data.load_data` for supported data formats and broadcasting
235 and reshaping rules.
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.
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[:]})
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)
257 # --------------------------------------------------------------------------
258 # (Abstract) class methods.
259 # --------------------------------------------------------------------------
261 @classmethod
262 @abstractmethod
263 def _get_bilinear_terms_allowed(cls):
264 """Report whether the subclass may have bilinear terms."""
265 pass
267 @classmethod
268 @abstractmethod
269 def _get_parameters_allowed(cls):
270 """Report whether the subclass may depend on parameters."""
271 pass
273 @classmethod
274 @abstractmethod
275 def _get_basetype(cls):
276 """Return first non-abstract :class:`Expression` subclass this bases on.
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
284 @classmethod
285 @abstractmethod
286 def _get_typecode(cls):
287 """Return the CVXOPT typecode to use with coefficient matrices.
289 Either ``"z"`` for complex or ``"d"`` for real.
291 See also :meth:`_get_basetype`.
292 """
293 pass
295 @classmethod
296 def _common_basetype(cls, other, reverse=False):
297 """Return the common basetype of two expressions.
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()
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__))
322 # --------------------------------------------------------------------------
323 # Internal use properties and methods.
324 # --------------------------------------------------------------------------
326 @property
327 def _bilinear_terms_allowed(self):
328 """Shorthand for :meth:`_get_bilinear_terms_allowed`."""
329 return self._get_bilinear_terms_allowed()
331 @property
332 def _parameters_allowed(self):
333 """Shorthand for :meth:`_get_parameters_allowed`."""
334 return self._get_parameters_allowed()
336 @property
337 def _basetype(self):
338 """Shorthand for :meth:`_get_basetype`."""
339 return self._get_basetype()
341 @property
342 def _typecode(self):
343 """Shorthand for :meth:`_get_typecode`."""
344 return self._get_typecode()
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]
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}
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}
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()}
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}
376 def _reglyphed_string(self, mutable_name_map):
377 """The symbolic string with mutable names replaced."""
378 string = self.string
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]
387 return string
389 # --------------------------------------------------------------------------
390 # Abstract method implementations and overridings for Expression.
391 # --------------------------------------------------------------------------
393 def _get_value(self):
394 # Create a copy of the constant term.
395 value = self._constant_coef[:]
397 # Add value of the linear part.
398 for mtb, coef in self._linear_coefs.items():
399 summand = coef * mtb._get_internal_value()
401 if type(value) is type(summand):
402 value += summand
403 else:
404 # Exactly one of the matrices is sparse.
405 value = value + summand
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)[:]
413 if type(value) is type(summand):
414 value += summand
415 else:
416 # Exactly one of the matrices is sparse.
417 value = value + summand
419 # Resize the value to the proper shape.
420 value.size = self.shape
422 return value
424 def _get_shape(self):
425 return self._shape
427 @cached_unary_operator
428 def _get_mutables(self):
429 return frozenset(mtb for mtbs in self._coefs for mtb in mtbs)
431 def _is_convex(self):
432 return not self._bilinear_coefs
434 def _is_concave(self):
435 return not self._bilinear_coefs
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]
442 name_map = {old.name: new.string for old, new in mapping.items()}
443 string = self._reglyphed_string(name_map)
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}
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
492 return self._basetype(string, self._shape, coefs)
494 def _freeze_mutables(self, freeze):
495 string = self._reglyphed_string(
496 {mtb.name: glyphs.frozen(mtb.name) for mtb in freeze})
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
517 return self._basetype(string, self._shape, coefs)
519 # --------------------------------------------------------------------------
520 # Python special method implementations.
521 # --------------------------------------------------------------------------
523 def __len__(self):
524 # Faster version that overrides Expression.__len__.
525 return self._shape[0] * self._shape[1]
527 def __getitem__(self, index):
528 def slice2range(s, length):
529 """Transform a :class:`slice` to a :class:`range`."""
530 assert isinstance(s, slice)
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.
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
550 # Clamp out-of-bound indices.
551 ra = min(max(0, ra), length - 1)
552 rb = min(max(-1, rb), length)
554 r = range(ra, rb, ss)
556 if not r:
557 raise IndexError("Empty slice.")
559 return r
561 def range2slice(r, length):
562 """Transform a :class:`range` to a :class:`slice`, if possible.
564 :raises ValueError: If the input cannot be expressed as a slice.
565 """
566 assert isinstance(r, range)
568 if not r:
569 raise IndexError("Empty range.")
571 ra = r.start
572 rb = r.stop
573 rs = r.step
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.")
585 if rb == -1:
586 rb = None
588 return slice(ra, rb, rs)
590 def list2slice(l, length):
591 """Transform a :class:`list` to a :class:`slice`, if possible.
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)
598 def slice2str(s):
599 """Return the short string that produced a :class:`slice`."""
600 assert isinstance(s, slice)
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))
613 def list2str(l, length):
614 """Return a short string represnetation of a :class:`list`."""
615 assert isinstance(l, list)
617 # Extract integers wrapped in a list.
618 if len(l) == 1:
619 return str(l[0])
621 # Try to convert the list to a slice.
622 try:
623 l = list2slice(l, length)
624 except (ValueError, RuntimeError):
625 pass
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)
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
643 m, n = self._shape
644 indexStr = None
645 isIntList = False
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.")
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
661 I, J = index[i], index[j]
663 try:
664 I = load_dense_data(I, typecode="i", alwaysCopy=False)[0]
665 J = load_dense_data(J, typecode="i", alwaysCopy=False)[0]
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.")
671 if len(I) != len(J):
672 raise TypeError("The objects do not have the same length.")
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
681 # Represent the selection as a global index list.
682 index = [[i + j*m for i, j in zip(I, J)]]
684 # Use a special index string.
685 indexStr = glyphs.size(list2str(I, n), list2str(J, m))
687 # Don't invoke load_dense_data on the list.
688 isIntList = True
689 else: # Global indexing.
690 index = [index]
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.")
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]
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)))
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
720 index[axis] = selection
722 # Build index string, retrieve new shape, finalize index.
723 if len(index) == 1: # Global indexing.
724 index = index[0]
726 if isinstance(index, slice):
727 shape = len(slice2range(index, len(self)))
728 else:
729 shape = len(index)
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))
738 # Convert to a global index list.
739 RC, shape = [], []
740 for axis, selection in enumerate(index):
741 k = self._shape[axis]
743 if isinstance(selection, slice):
744 # Turn the slice into an iterable range.
745 selection = slice2range(selection, k)
747 # All indices from a slice are nonnegative.
748 assert all(i >= 0 for i in selection)
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]
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))
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))
764 RC.append(selection)
765 shape.append(len(selection))
767 rows, cols = RC
768 index = [i + j*m for j in cols for i in rows]
770 # Finalize the string.
771 string = glyphs.slice(self.string, indexStr)
773 # Retrieve new coefficients and constant term.
774 coefs = {mtbs: coef[index, :] for mtbs, coef in self._coefs.items()}
776 return self._basetype(string, shape, coefs)
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)
784 string = glyphs.clever_add(self.string, other.string)
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)
793 return self._common_basetype(other)(string, self._shape, coefs)
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)
801 return other.__add__(self)
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)
809 string = glyphs.clever_sub(self.string, other.string)
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)
818 return self._common_basetype(other)(string, self._shape, coefs)
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)
826 return other.__sub__(self)
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()}
833 return self._basetype(string, self._shape, coefs)
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)
841 string = glyphs.clever_mul(self.string, other.string)
843 m, n = self._shape
844 p, q = other._shape
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
857 shape = other._shape if m*n == 1 else self._shape
859 # Work around CVXOPT not broadcasting sparse scalars.
860 if lhs.size == (1, 1):
861 lhs = lhs[0]
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
879 assert n == p
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
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))
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)
910 self._update_coefs(coefs, lmtbs, factor*lcoef)
911 else:
912 # Linear by linear.
913 assert len(lmtbs) == 1 and len(rmtbs) == 1
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
923 A = cvxopt_K(m, n)*lcoef
924 B = rcoef
926 A = A.T
927 B = B.T
928 A.size = m*n*a, 1
929 B.size = p*q*b, 1
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
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
943 coef = C
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))
951 self._update_coefs(coefs, lmtbs + rmtbs, coef)
953 return self._common_basetype(other)(string, (m, q), coefs)
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)
961 return other.__mul__(self)
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)
969 result = self.vec.H * other.vec
970 result._symbStr = glyphs.clever_dotp(
971 self.string, other.string, self.complex, self.scalar)
972 return result
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)
980 return other.__or__(self)
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)
988 string = glyphs.hadamard(self.string, other.string)
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)
999 return self._common_basetype(other)(string, self._shape, coefs)
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)
1007 return other.__xor__(self)
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)
1015 cls = self._common_basetype(other)
1016 tc = cls._get_typecode()
1018 string = glyphs.kron(self.string, other.string)
1020 m, n = self._shape
1021 p, q = other._shape
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)
1028 def _kron(A, B):
1029 if isinstance(A, cvxopt.spmatrix) or isinstance(B, cvxopt.spmatrix):
1030 try:
1031 import scipy
1032 except ModuleNotFoundError:
1033 A_B = load_data(
1034 numpy.kron(cvx2np(A), cvx2np(B)), typecode=tc
1035 )[0]
1036 else:
1037 from .data import cvx2csc
1038 A_B = load_data(
1039 scipy.sparse.kron(cvx2csc(A), cvx2csc(B), format='csc'),
1040 typecode=tc, sparse=True
1041 )[0]
1042 else:
1043 A_B = load_data(
1044 numpy.kron(cvx2np(A), cvx2np(B)), typecode=tc
1045 )[0]
1046 Kab = cvxopt_K(A.size[1], B.size[1])
1047 return In_Kqm_Ip*A_B*Kab
1049 coefs = {}
1050 for (lmtbs, lcoef), (rmtbs, rcoef) in product(
1051 self._coefs.items(), other._coefs.items()):
1052 if len(lmtbs) + len(rmtbs) <= 2:
1053 if len(lmtbs) == 1 and len(rmtbs) == 1 and lmtbs[0] is rmtbs[0]:
1054 raise TypeError("The product of {} and {} is quadratic "
1055 "in {} but a biaffine result is required here."
1056 .format(self.string, other.string, lmtbs[0].name))
1058 self._update_coefs(coefs, lmtbs + rmtbs, _kron(lcoef, rcoef))
1059 else:
1060 raise NotImplementedError("Product not supported: "
1061 "One operand is biaffine and the other nonconstant.")
1063 return cls(string, (m*p, n*q), coefs)
1065 @convert_operands()
1066 @refine_operands(stop_at_affine=True)
1067 def __rmatmul__(self, other):
1068 if not isinstance(other, BiaffineExpression):
1069 return Expression.__rmatmul__(self, other)
1071 return other.__matmul__(self)
1073 @deprecated("2.2", "Use infix @ instead.")
1074 def kron(self, other):
1075 """Denote the Kronecker product with another expression on the right."""
1076 return self.__matmul__(other)
1078 @deprecated("2.2", "Reverse operands and use infix @ instead.")
1079 def leftkron(self, other):
1080 """Denote the Kronecker product with another expression on the left."""
1081 return self.__rmatmul__(other)
1083 @convert_operands(scalarRHS=True)
1084 @refine_operands(stop_at_affine=True)
1085 def __truediv__(self, other):
1086 if not isinstance(other, BiaffineExpression):
1087 return Expression.__truediv__(self, other)
1089 if not other.constant:
1090 raise TypeError(
1091 "You may only divide {} by a constant.".format(self.string))
1093 if other.is0:
1094 raise ZeroDivisionError(
1095 "Tried to divide {} by zero.".format(self.string))
1097 divisor = other._constant_coef[0]
1099 string = glyphs.div(self.string, other.string)
1100 coefs = {mtbs: coef / divisor for mtbs, coef in self._coefs.items()}
1102 return self._common_basetype(other)(string, self._shape, coefs)
1104 @convert_operands(scalarLHS=True)
1105 @refine_operands(stop_at_affine=True)
1106 def __rtruediv__(self, other):
1107 if not isinstance(other, BiaffineExpression):
1108 return Expression.__rtruediv__(self, other)
1110 return other.__truediv__(self)
1112 @convert_operands(horiCat=True)
1113 @refine_operands(stop_at_affine=True)
1114 def __and__(self, other):
1115 if not isinstance(other, BiaffineExpression):
1116 return Expression.__and__(self, other)
1118 string = glyphs.matrix_cat(self.string, other.string, horizontal=True)
1119 shape = (self._shape[0], self._shape[1] + other._shape[1])
1121 coefs = {}
1122 for mtbs in set(self._coefs.keys()).union(other._coefs.keys()):
1123 l = self._coefs[mtbs] if mtbs in self._coefs else cvxopt.spmatrix(
1124 [], [], [], (len(self), other._coefs[mtbs].size[1]))
1125 r = other._coefs[mtbs] if mtbs in other._coefs else cvxopt.spmatrix(
1126 [], [], [], (len(other), self._coefs[mtbs].size[1]))
1128 coefs[mtbs] = cvxopt_vcat([l, r])
1130 return self._common_basetype(other)(string, shape, coefs)
1132 @convert_operands(horiCat=True)
1133 @refine_operands(stop_at_affine=True)
1134 def __rand__(self, other):
1135 if not isinstance(other, BiaffineExpression):
1136 return Expression.__rand__(self, other)
1138 return other.__and__(self)
1140 @convert_operands(vertCat=True)
1141 @refine_operands(stop_at_affine=True)
1142 def __floordiv__(self, other):
1143 def interleave_columns(upper, lower, upperRows, lowerRows, cols):
1144 p, q = upperRows, lowerRows
1145 return [column for columnPairs in [
1146 (upper[j*p:j*p+p, :], lower[j*q:j*q+q, :]) for j in range(cols)]
1147 for column in columnPairs]
1149 if not isinstance(other, BiaffineExpression):
1150 return Expression.__floordiv__(self, other)
1152 string = glyphs.matrix_cat(self.string, other.string, horizontal=False)
1154 p, q, c = self._shape[0], other._shape[0], self._shape[1]
1155 shape = (p + q, c)
1157 coefs = {}
1158 for mtbs in set(self._coefs.keys()).union(other._coefs.keys()):
1159 u = self._coefs[mtbs] if mtbs in self._coefs else cvxopt.spmatrix(
1160 [], [], [], (len(self), other._coefs[mtbs].size[1]))
1161 l = other._coefs[mtbs] if mtbs in other._coefs else cvxopt.spmatrix(
1162 [], [], [], (len(other), self._coefs[mtbs].size[1]))
1164 coefs[mtbs] = cvxopt_vcat(interleave_columns(u, l, p, q, c))
1166 return self._common_basetype(other)(string, shape, coefs)
1168 @convert_operands(vertCat=True)
1169 @refine_operands(stop_at_affine=True)
1170 def __rfloordiv__(self, other):
1171 if not isinstance(other, BiaffineExpression):
1172 return Expression.__rfloordiv__(self, other)
1174 return other.__floordiv__(self)
1176 @convert_operands(scalarRHS=True)
1177 @refine_operands() # Refine both sides to real if possible.
1178 def __pow__(self, other):
1179 from .exp_powtrace import PowerTrace
1181 if not self.scalar:
1182 raise TypeError("May only exponentiate a scalar expression.")
1184 if not other.constant:
1185 raise TypeError("The exponent must be constant.")
1187 if other.complex:
1188 raise TypeError("The exponent must be real-valued.")
1190 exponent = other.value
1192 if exponent == 2:
1193 return (self | self) # Works for complex base.
1194 else:
1195 return PowerTrace(self, exponent) # Errors on complex base.
1197 # --------------------------------------------------------------------------
1198 # Properties and functions that describe the expression.
1199 # --------------------------------------------------------------------------
1201 @cached_property
1202 def hermitian(self): # noqa (D402 thinks this includes a signature)
1203 """Whether the expression is a hermitian (or symmetric) matrix.
1205 Uses :data:`~picos.settings.RELATIVE_HERMITIANNESS_TOLERANCE`.
1207 If PICOS rejects your near-hermitian (near-symmetric) expression as not
1208 hermitian (not symmetric), you can use :meth:`hermitianized` to correct
1209 larger numeric errors or the effects of noisy data.
1210 """
1211 return self.equals(
1212 self.H, relTol=settings.RELATIVE_HERMITIANNESS_TOLERANCE)
1214 @property
1215 def is0(self):
1216 """Whether this is a constant scalar, vector or matrix of all zeros."""
1217 return not self._coefs
1219 @cached_property
1220 def is1(self):
1221 """Whether this is a constant scalar or vector of all ones."""
1222 # Must be a scalar or vector.
1223 if self._shape[0] != 1 and self._shape[1] != 1:
1224 return False
1226 # Must be constant.
1227 if not self.constant:
1228 return False
1230 # Constant term must be all ones.
1231 return not self._constant_coef - 1
1233 @cached_property
1234 def isI(self):
1235 """Whether this is a constant identity matrix."""
1236 m, n = self._shape
1238 # Must be a square matrix.
1239 if m != n:
1240 return False
1242 # Must be constant.
1243 if not self.constant:
1244 return False
1246 # Constant term must be the identity.
1247 return not self._constant_coef - cvxopt.spdiag([1]*m)[:]
1249 @cached_property
1250 def complex(self):
1251 """Whether the expression can be complex-valued."""
1252 # NOTE: The typecode check works around a bug in CVXOPT prior to 1.2.4.
1253 return any(coef.typecode == "z" and coef.imag()
1254 for coef in self._coefs.values())
1256 @property
1257 def isreal(self):
1258 """Whether the expression is always real-valued."""
1259 return not self.complex
1261 @convert_operands()
1262 def equals(self, other, absTol=None, relTol=None):
1263 """Check mathematical equality with another (bi)affine expression.
1265 The precise type of both (bi)affine expressions may differ. In
1266 particular, a :class:`~.exp_affine.ComplexAffineExpression` with real
1267 coefficients and constant term can be equal to an
1268 :class:`~.exp_affine.AffineExpression`.
1270 If the operand is not already a PICOS expression, an attempt is made to
1271 load it as a constant affine expression. In this case, no reshaping or
1272 broadcasting is used to bring the constant term to the same shape as
1273 this expression. In particular,
1275 - ``0`` refers to a scalar zero (see also :meth:`is0`),
1276 - lists and tuples are treated as column vectors and
1277 - algebraic strings must specify a shape (see
1278 :func:`~.data.load_data`).
1280 :param other: Another PICOS expression or a constant numeric data value
1281 supported by :func:`~.data.load_data`.
1282 :param absTol: As long as all absolute differences between scalar
1283 entries of the coefficient matrices and the constant terms being
1284 compared does not exceed this bound, consider the expressions equal.
1285 :param relTol: As long as all absolute differences between scalar
1286 entries of the coefficient matrices and the constant terms being
1287 compared divided by the maximum absolute value found in either term
1288 does not exceed this bound, consider the expressions equal.
1290 :Example:
1292 >>> from picos import Constant
1293 >>> A = Constant("A", 0, (5,5))
1294 >>> repr(A)
1295 '<5×5 Real Constant: A>'
1296 >>> A.is0
1297 True
1298 >>> A.equals(0)
1299 False
1300 >>> A.equals("|0|(5,5)")
1301 True
1302 >>> repr(A*1j)
1303 '<5×5 Complex Constant: A·1j>'
1304 >>> A.equals(A*1j)
1305 True
1306 """
1307 if self is other:
1308 return True
1310 if not isinstance(other, BiaffineExpression):
1311 return False
1313 if self._shape != other._shape:
1314 return False
1316 assert not any((mtbs[1], mtbs[0]) in other._bilinear_coefs
1317 for mtbs in self._bilinear_coefs), \
1318 "{} and {} store bilinear terms in a different order." \
1319 .format(self.string, other.string)
1321 for mtbs in other._coefs:
1322 if mtbs not in self._coefs:
1323 return False
1325 for mtbs, coef in self._coefs.items():
1326 if mtbs not in other._coefs:
1327 return False
1329 if not cvxopt_equals(coef, other._coefs[mtbs], absTol, relTol):
1330 return False
1332 return True
1334 # --------------------------------------------------------------------------
1335 # Methods and properties that return modified copies.
1336 # --------------------------------------------------------------------------
1338 def renamed(self, string):
1339 """Return the expression with a modified string description."""
1340 return self._basetype(string, self._shape, self._coefs)
1342 def reshaped(self, shape, order="F"):
1343 r"""Return the expression reshaped in the given order.
1345 The default indexing order is column-major. Given an :math:`m \times n`
1346 matrix, reshaping in the default order is a constant time operation
1347 while reshaping in row-major order requires :math:`O(mn)` time. However,
1348 the latter allows you to be consistent with NumPy, which uses C-order (a
1349 generalization of row-major) by default.
1351 :param str order:
1352 The indexing order to use when reshaping. Must be either ``"F"`` for
1353 Fortran-order (column-major) or ``"C"`` for C-order (row-major).
1355 :Example:
1357 >>> from picos import Constant
1358 >>> C = Constant("C", range(6), (2, 3))
1359 >>> print(C)
1360 [ 0.00e+00 2.00e+00 4.00e+00]
1361 [ 1.00e+00 3.00e+00 5.00e+00]
1362 >>> print(C.reshaped((3, 2)))
1363 [ 0.00e+00 3.00e+00]
1364 [ 1.00e+00 4.00e+00]
1365 [ 2.00e+00 5.00e+00]
1366 >>> print(C.reshaped((3, 2), order="C"))
1367 [ 0.00e+00 2.00e+00]
1368 [ 4.00e+00 1.00e+00]
1369 [ 3.00e+00 5.00e+00]
1370 """
1371 if order not in "FC":
1372 raise ValueError("Order must be given as 'F' or 'C'.")
1374 shape = load_shape(shape, wildcards=True)
1376 if shape == self._shape:
1377 return self
1378 elif shape == (None, None):
1379 return self
1381 length = len(self)
1383 if shape[0] is None:
1384 shape = (length // shape[1], shape[1])
1385 elif shape[1] is None:
1386 shape = (shape[0], length // shape[0])
1388 if shape[0]*shape[1] != length:
1389 raise ValueError("Can only reshape to a matrix of same size.")
1391 if order == "F":
1392 coefs = self._coefs
1393 reshaped_glyph = glyphs.reshaped
1394 else:
1395 m, n = self._shape
1396 p, q = shape
1398 K_old = cvxopt_K(m, n, self._typecode)
1399 K_new = cvxopt_K(q, p, self._typecode)
1400 R = K_new * K_old
1402 coefs = {mtbs: R * coef for mtbs, coef in self._coefs.items()}
1403 reshaped_glyph = glyphs.reshaprm
1405 string = reshaped_glyph(self.string, glyphs.shape(shape))
1406 return self._basetype(string, shape, coefs)
1408 def broadcasted(self, shape):
1409 """Return the expression broadcasted to the given shape.
1411 :Example:
1413 >>> from picos import Constant
1414 >>> C = Constant("C", range(6), (2, 3))
1415 >>> print(C)
1416 [ 0.00e+00 2.00e+00 4.00e+00]
1417 [ 1.00e+00 3.00e+00 5.00e+00]
1418 >>> print(C.broadcasted((6, 6)))
1419 [ 0.00e+00 2.00e+00 4.00e+00 0.00e+00 2.00e+00 4.00e+00]
1420 [ 1.00e+00 3.00e+00 5.00e+00 1.00e+00 3.00e+00 5.00e+00]
1421 [ 0.00e+00 2.00e+00 4.00e+00 0.00e+00 2.00e+00 4.00e+00]
1422 [ 1.00e+00 3.00e+00 5.00e+00 1.00e+00 3.00e+00 5.00e+00]
1423 [ 0.00e+00 2.00e+00 4.00e+00 0.00e+00 2.00e+00 4.00e+00]
1424 [ 1.00e+00 3.00e+00 5.00e+00 1.00e+00 3.00e+00 5.00e+00]
1425 """
1426 shape = load_shape(shape, wildcards=True)
1427 shape = blend_shapes(shape, self._shape)
1429 if shape == self._shape:
1430 return self
1432 vdup = shape[0] // self._shape[0]
1433 hdup = shape[1] // self._shape[1]
1435 if (self._shape[0] * vdup, self._shape[1] * hdup) != shape:
1436 raise ValueError("Cannot broadcast from shape {} to {}."
1437 .format(glyphs.shape(self._shape), glyphs.shape(shape)))
1439 if self._shape == (1, 1):
1440 string = glyphs.matrix(self.string)
1441 return (self * cvxopt.matrix(1.0, shape)).renamed(string)
1443 string = glyphs.bcasted(self.string, glyphs.shape(shape))
1444 return (cvxopt.matrix(1.0, (vdup, hdup)) @ self).renamed(string)
1446 def reshaped_or_broadcasted(self, shape):
1447 """Return the expression :meth:`reshaped` or :meth:`broadcasted`.
1449 Unlike with :meth:`reshaped` and :meth:`broadcasted`, the target shape
1450 may not contain a wildcard character.
1452 If the wildcard-free target shape has the same number of elements as
1453 the current shape, then this is the same as :meth:`reshaped`, otherwise
1454 it is the same as :meth:`broadcasted`.
1455 """
1456 shape = load_shape(shape, wildcards=False)
1458 try:
1459 if shape[0]*shape[1] == len(self):
1460 return self.reshaped(shape)
1461 else:
1462 return self.broadcasted(shape)
1463 except ValueError:
1464 raise ValueError(
1465 "Cannot reshape or broadcast from shape {} to {}.".format(
1466 glyphs.shape(self._shape), glyphs.shape(shape))) from None
1468 @cached_property
1469 def opreal(self):
1470 r"""The real part of a Hermitian operator.
1472 For a square (complex) affine expression :math:`A`, this is
1473 :math:`\frac{1}{2}(A + A^H)`.
1475 :Example:
1477 >>> from picos import ComplexVariable
1478 >>> ComplexVariable("X", (4, 4)).opreal.hermitian
1479 True
1480 """
1481 if not self.square:
1482 raise TypeError("Cannot take opreal of non-square {}.".format(self))
1484 return (self + self.H)/2
1486 @cached_property
1487 def opimag(self):
1488 r"""The imaginary part of a Hermitian operator.
1490 For a square (complex) affine expression :math:`A`, this is
1491 :math:`\frac{1}{2i}(A - A^H)`.
1493 :Example:
1495 >>> from picos import ComplexVariable
1496 >>> ComplexVariable("X", (4, 4)).opimag.hermitian
1497 True
1498 """
1499 if not self.square:
1500 raise TypeError("Cannot take opimag of non-square {}.".format(self))
1502 return (self - self.H)/(2j)
1504 @property
1505 def hermitianized(self):
1506 r"""The expression projected onto the subspace of hermitian matrices.
1508 If the expression is not complex, then this is a projection onto the
1509 subspace of symmetric matrices.
1511 This is the same as :meth:`opreal`.
1512 """
1513 return self.opreal
1515 @cached_property
1516 def real(self):
1517 """Real part of the expression."""
1518 return self._basetype(glyphs.real(self.string), self._shape,
1519 {mtbs: coef.real() for mtbs, coef in self._coefs.items()})
1521 @cached_property
1522 def imag(self):
1523 """Imaginary part of the expression."""
1524 return self._basetype(glyphs.imag(self.string), self._shape,
1525 {mtbs: coef.imag() for mtbs, coef in self._coefs.items()})
1527 @cached_property
1528 def bilin(self):
1529 """Pure bilinear part of the expression."""
1530 return self._basetype(
1531 glyphs.blinpart(self._symbStr), self._shape, self._bilinear_coefs)
1533 @cached_property
1534 def lin(self):
1535 """Linear part of the expression."""
1536 return self._basetype(
1537 glyphs.linpart(self._symbStr), self._shape, self._linear_coefs)
1539 @cached_property
1540 def noncst(self):
1541 """Nonconstant part of the expression."""
1542 coefs = {mtbs: coefs for mtbs, coefs in self._coefs.items() if mtbs}
1544 return self._basetype(
1545 glyphs.ncstpart(self._symbStr), self._shape, coefs)
1547 @cached_property
1548 def cst(self):
1549 """Constant part of the expression."""
1550 coefs = {(): self._coefs[()]} if () in self._coefs else {}
1552 return self._basetype(glyphs.cstpart(self._symbStr), self._shape, coefs)
1554 @cached_selfinverse_property
1555 def T(self):
1556 """Matrix transpose."""
1557 if len(self) == 1:
1558 return self
1560 m, n = self._shape
1561 K = cvxopt_K(m, n, self._typecode)
1563 string = glyphs.transp(self.string)
1564 shape = (self._shape[1], self._shape[0])
1565 coefs = {mtbs: K * coef for mtbs, coef in self._coefs.items()}
1567 return self._basetype(string, shape, coefs)
1569 @cached_selfinverse_property
1570 def conj(self):
1571 """Complex conjugate."""
1572 string = glyphs.conj(self.string)
1573 coefs = {mtbs: coef.H.T for mtbs, coef in self._coefs.items()}
1575 return self._basetype(string, self._shape, coefs)
1577 @cached_selfinverse_property
1578 def H(self):
1579 """Conjugate (or Hermitian) transpose."""
1580 return self.T.conj.renamed(glyphs.htransp(self._symbStr))
1582 def _square_equal_subsystem_dims(self, diagLen):
1583 """Support :func:`partial_trace` and :func:`partial_transpose`."""
1584 m, n = self._shape
1585 k = math.log(m, diagLen)
1587 if m != n or int(k) != k:
1588 raise TypeError("The expression has shape {} so it cannot be "
1589 "decomposed into subsystems of shape {}.".format(
1590 glyphs.shape(self._shape), glyphs.shape((diagLen,)*2)))
1592 return ((diagLen,)*2,)*int(k)
1594 def partial_transpose(self, subsystems, dimensions=2):
1595 r"""Return the expression with selected subsystems transposed.
1597 If the expression can be written as
1598 :math:`A_0 \otimes \cdots \otimes A_{n-1}` for matrices
1599 :math:`A_0, \ldots, A_{n-1}` with shapes given in ``dimensions``, then
1600 this returns :math:`B_0 \otimes \cdots \otimes B_{n-1}` with
1601 :math:`B_i = A_i^T`, if ``i in subsystems`` (with :math:`i = -1` read as
1602 :math:`n-1`), and :math:`B_i = A_i`, otherwise.
1604 :param subsystems: A collection of or a single subystem number, indexed
1605 from zero, corresponding to subsystems that shall be transposed.
1606 The value :math:`-1` refers to the last subsystem.
1607 :type subsystems: int or tuple or list
1609 :param dimensions: Either an integer :math:`d` so that the subsystems
1610 are assumed to be all of shape :math:`d \times d`, or a sequence of
1611 subsystem shapes where an integer :math:`d` within the sequence is
1612 read as :math:`d \times d`. In any case, the elementwise product
1613 over all subsystem shapes must equal the expression's shape.
1614 :type dimensions: int or tuple or list
1616 :raises TypeError: If the subsystems do not match the expression.
1617 :raises IndexError: If the subsystem selection is invalid.
1619 :Example:
1621 >>> from picos import Constant
1622 >>> A = Constant("A", range(16), (4, 4))
1623 >>> print(A) #doctest: +NORMALIZE_WHITESPACE
1624 [ 0.00e+00 4.00e+00 8.00e+00 1.20e+01]
1625 [ 1.00e+00 5.00e+00 9.00e+00 1.30e+01]
1626 [ 2.00e+00 6.00e+00 1.00e+01 1.40e+01]
1627 [ 3.00e+00 7.00e+00 1.10e+01 1.50e+01]
1628 >>> A0 = A.partial_transpose(0); A0
1629 <4×4 Real Constant: A.{[2×2]ᵀ⊗[2×2]}>
1630 >>> print(A0) #doctest: +NORMALIZE_WHITESPACE
1631 [ 0.00e+00 4.00e+00 2.00e+00 6.00e+00]
1632 [ 1.00e+00 5.00e+00 3.00e+00 7.00e+00]
1633 [ 8.00e+00 1.20e+01 1.00e+01 1.40e+01]
1634 [ 9.00e+00 1.30e+01 1.10e+01 1.50e+01]
1635 >>> A1 = A.partial_transpose(1); A1
1636 <4×4 Real Constant: A.{[2×2]⊗[2×2]ᵀ}>
1637 >>> print(A1) #doctest: +NORMALIZE_WHITESPACE
1638 [ 0.00e+00 1.00e+00 8.00e+00 9.00e+00]
1639 [ 4.00e+00 5.00e+00 1.20e+01 1.30e+01]
1640 [ 2.00e+00 3.00e+00 1.00e+01 1.10e+01]
1641 [ 6.00e+00 7.00e+00 1.40e+01 1.50e+01]
1642 """
1643 m, n = self._shape
1645 if isinstance(dimensions, int):
1646 dimensions = self._square_equal_subsystem_dims(dimensions)
1647 else:
1648 dimensions = [
1649 (d, d) if isinstance(d, int) else d for d in dimensions]
1651 if reduce(
1652 lambda x, y: (x[0]*y[0], x[1]*y[1]), dimensions) != self._shape:
1653 raise TypeError("Subsystem dimensions do not match expression.")
1655 if isinstance(subsystems, int):
1656 subsystems = (subsystems,)
1657 elif not subsystems:
1658 return self
1660 numSys = len(dimensions)
1661 subsystems = set(numSys - 1 if sys == -1 else sys for sys in subsystems)
1663 for sys in subsystems:
1664 if not isinstance(sys, int):
1665 raise IndexError("Subsystem indices must be integer, not {}."
1666 .format(type(sys).__name__))
1667 elif sys < 0:
1668 raise IndexError("Subsystem indices must be nonnegative.")
1669 elif sys >= numSys:
1670 raise IndexError(
1671 "Subsystem index {} out of range for {} systems total."
1672 .format(sys, numSys))
1674 # If all subsystems are transposed, this is regular transposition.
1675 if len(subsystems) == numSys:
1676 return self.T
1678 # Prepare sparse K such that K·vec(A) = vec(partial_transpose(A)).
1679 d = m * n
1680 V = [1]*d
1681 I = range(d)
1682 J = cvxopt.matrix(I)
1683 T = cvxopt.matrix(0, J.size)
1684 obh, obw = 1, 1
1685 sysStrings = None
1687 # Apply transpositions starting with the rightmost Kronecker factor.
1688 for sys in range(numSys - 1, -1, -1):
1689 # Shape of current system.
1690 p, q = dimensions[sys]
1691 sysString = glyphs.matrix(glyphs.shape((p, q)))
1693 # Height/width of "inner" blocks being moved, initially scalars.
1694 ibh, ibw = obh, obw
1696 # Heigh/width of "outer" blocks whose relative position is
1697 # maintained but that are subject to transposition independently.
1698 # In the last iteration this is the shape of the resulting matrix.
1699 obh *= p
1700 obw *= q
1702 # Only transpose selected subsystems.
1703 if sys not in subsystems:
1704 sysStrings = glyphs.kron(sysString, sysStrings) \
1705 if sysStrings else sysString
1706 continue
1707 else:
1708 sysStrings = glyphs.kron(glyphs.transp(sysString), sysStrings) \
1709 if sysStrings else glyphs.transp(sysString)
1711 # Shape of outer blocks after transposition.
1712 obhT, obwT = obw // ibw * ibh, obh // ibh * ibw
1714 # Shape of full matrix after transposition.
1715 mT, nT = m // obh * obhT, n // obw * obwT
1717 for vi in I:
1718 # Full matrix column and row.
1719 c, r = divmod(vi, m)
1721 # Outer block vertical index and row within outer block,
1722 # outer block horizontal index and column within outer block.
1723 obi, obr = divmod(r, obh)
1724 obj, obc = divmod(c, obw)
1726 # Inner block vertical index and row within inner block,
1727 # inner block horizontal index and column within inner block.
1728 ibi, ibr = divmod(obr, ibh)
1729 ibj, ibc = divmod(obc, ibw)
1731 # (1) ibi*ibw + ibc is column within the transposed outer block;
1732 # adding obj*obwT yields the column in the transposed matrix.
1733 # (2) ibj*ibh + ibr is row within the transposed outer block;
1734 # adding obi*obhT yields the row in the transposed matrix.
1735 # (3) tvi is index within the vectorized transposed matrix.
1736 tvi = (obj*obwT + ibi*ibw + ibc)*mT \
1737 + (obi*obhT + ibj*ibh + ibr)
1739 # Prepare the transposition.
1740 T[tvi] = J[vi]
1742 # Apply the transposition.
1743 J, T = T, J
1744 m, n, obh, obw = mT, nT, obhT, obwT
1746 # Finalize the partial commutation matrix.
1747 K = cvxopt.spmatrix(V, I, J, (d, d), self._typecode)
1749 string = glyphs.ptransp_(self.string, sysStrings)
1750 shape = (m, n)
1751 coefs = {mtbs: K * coef for mtbs, coef in self._coefs.items()}
1753 return self._basetype(string, shape, coefs)
1755 @cached_property
1756 def T0(self):
1757 r"""Expression with the first :math:`2 \times 2` subsystem transposed.
1759 Only available for a :math:`2^k \times 2^k` matrix with all subsystems
1760 of shape :math:`2 \times 2`. Use :meth:`partial_transpose` otherwise.
1761 """
1762 return self.partial_transpose(subsystems=0)
1764 @cached_property
1765 def T1(self):
1766 r"""Expression with the second :math:`2 \times 2` subsystem transposed.
1768 Only available for a :math:`2^k \times 2^k` matrix with all subsystems
1769 of shape :math:`2 \times 2`. Use :meth:`partial_transpose` otherwise.
1770 """
1771 return self.partial_transpose(subsystems=1)
1773 @cached_property
1774 def T2(self):
1775 r"""Expression with the third :math:`2 \times 2` subsystem transposed.
1777 Only available for a :math:`2^k \times 2^k` matrix with all subsystems
1778 of shape :math:`2 \times 2`. Use :meth:`partial_transpose` otherwise.
1779 """
1780 return self.partial_transpose(subsystems=2)
1782 @cached_property
1783 def T3(self):
1784 r"""Expression with the fourth :math:`2 \times 2` subsystem transposed.
1786 Only available for a :math:`2^k \times 2^k` matrix with all subsystems
1787 of shape :math:`2 \times 2`. Use :meth:`partial_transpose` otherwise.
1788 """
1789 return self.partial_transpose(subsystems=3)
1791 @cached_property
1792 def Tl(self):
1793 r"""Expression with the last :math:`2 \times 2` subsystem transposed.
1795 Only available for a :math:`2^k \times 2^k` matrix with all subsystems
1796 of shape :math:`2 \times 2`. Use :meth:`partial_transpose` otherwise.
1797 """
1798 return self.partial_transpose(subsystems=-1)
1800 @staticmethod
1801 def _reindex_F(indices, source, destination):
1802 """Convert indices between different tensor shapes in Fortran-order."""
1803 new = []
1804 offset = 0
1805 factor = 1
1807 for index, base in zip(indices, source):
1808 offset += factor*index
1809 factor *= base
1811 for base in destination:
1812 offset, remainder = divmod(offset, base)
1813 new.append(remainder)
1815 return tuple(new)
1817 @staticmethod
1818 def _reindex_C(indices, source, destination):
1819 """Convert indices between different tensor shapes in C-order."""
1820 new = []
1821 offset = 0
1822 factor = 1
1824 for index, base in zip(reversed(indices), reversed(source)):
1825 offset += factor*index
1826 factor *= base
1828 for base in reversed(destination):
1829 offset, remainder = divmod(offset, base)
1830 new.insert(0, remainder)
1832 return tuple(new)
1834 def reshuffled(self, permutation="ikjl", dimensions=None, order="C"):
1835 """Return the reshuffled or realigned expression.
1837 This operation works directly on matrices. However, it is equivalent to
1838 the following sequence of operations:
1840 1. The matrix is reshaped to a tensor with the given ``dimensions`` and
1841 according to ``order``.
1842 2. The tensor's axes are permuted according to ``permutation``.
1843 3. The tensor is reshaped back to the shape of the original matrix
1844 according to ``order``.
1846 For comparison, the following function applies the same operation to a
1847 2D NumPy :class:`~numpy:numpy.ndarray`:
1849 .. code::
1851 def reshuffle_numpy(matrix, permutation, dimensions, order):
1852 P = "{} -> {}".format("".join(sorted(permutation)), permutation)
1853 reshuffled = numpy.reshape(matrix, dimensions, order)
1854 reshuffled = numpy.einsum(P, reshuffled)
1855 return numpy.reshape(reshuffled, matrix.shape, order)
1857 :param permutation:
1858 A sequence of comparable elements with length equal to the number of
1859 tensor dimensions. The sequence is compared to its ordered version
1860 and the resulting permutation pattern is used to permute the tensor
1861 indices. For instance, the string ``"ikjl"`` is compared to its
1862 sorted version ``"ijkl"`` and denotes that the second and third axis
1863 should be swapped.
1864 :type permutation:
1865 str or tuple or list
1867 :param dimensions:
1868 If this is an integer sequence, then it defines the dimensions of
1869 the tensor. If this is :obj:`None`, then the tensor is assumed to be
1870 hypercubic and the number of dimensions is inferred from the
1871 ``permutation`` argument.
1872 :type dimensions:
1873 None or tuple or list
1875 :param str order:
1876 The indexing order to use for the virtual reshaping. Must be either
1877 ``"F"`` for Fortran-order (generalization of column-major) or
1878 ``"C"`` for C-order (generalization of row-major). Note that PICOS
1879 usually reshapes in Fortran-order while NumPy defaults to C-order.
1881 :Example:
1883 >>> from picos import Constant
1884 >>> A = Constant("A", range(16), (4, 4))
1885 >>> print(A) #doctest: +NORMALIZE_WHITESPACE
1886 [ 0.00e+00 4.00e+00 8.00e+00 1.20e+01]
1887 [ 1.00e+00 5.00e+00 9.00e+00 1.30e+01]
1888 [ 2.00e+00 6.00e+00 1.00e+01 1.40e+01]
1889 [ 3.00e+00 7.00e+00 1.10e+01 1.50e+01]
1890 >>> R = A.reshuffled(); R
1891 <4×4 Real Constant: shuffled(A,ikjl,C)>
1892 >>> print(R) #doctest: +NORMALIZE_WHITESPACE
1893 [ 0.00e+00 4.00e+00 1.00e+00 5.00e+00]
1894 [ 8.00e+00 1.20e+01 9.00e+00 1.30e+01]
1895 [ 2.00e+00 6.00e+00 3.00e+00 7.00e+00]
1896 [ 1.00e+01 1.40e+01 1.10e+01 1.50e+01]
1897 >>> A.reshuffled("ji").equals(A.T) # Regular transposition.
1898 True
1899 >>> A.reshuffled("3214").equals(A.T0) # Partial transposition (1).
1900 True
1901 >>> A.reshuffled("1432").equals(A.T1) # Partial transposition (2).
1902 True
1903 """
1904 m, n = self._shape
1905 mn = m*n
1907 # Load the permutation.
1908 ordered = sorted(permutation)
1909 P = dict(enumerate(ordered.index(element) for element in permutation))
1911 if len(set(P.values())) < len(P):
1912 raise ValueError("The sequence defining the permutation appears to "
1913 "contain duplicate elements.")
1915 assert not set(P.keys()).symmetric_difference(set(P.values()))
1917 numDims = len(P)
1919 # Load the dimensions.
1920 guessDimensions = dimensions is None
1922 if guessDimensions:
1923 dimensions = (int(mn**(1.0 / numDims)),)*numDims
1924 else:
1925 if len(dimensions) != numDims:
1926 raise ValueError("The number of indices does not match the "
1927 "number of dimensions.")
1929 if reduce(int.__mul__, dimensions, 1) != mn:
1930 raise TypeError("The {} matrix {} cannot be reshaped to a {} "
1931 "tensor.".format(glyphs.shape(self.shape), self.string,
1932 "hypercubic order {}".format(numDims) if guessDimensions
1933 else glyphs.size("", "").join(str(d) for d in dimensions)))
1935 # Load the indexing order.
1936 if order not in "FC":
1937 raise ValueError("Order must be given as 'F' or 'C'.")
1939 reindex = self._reindex_F if order == "F" else self._reindex_C
1941 # Nothing to do for the neutral permutation.
1942 if all(key == val for key, val in P.items()):
1943 return self
1945 # Create a sparse mtrix R such that R·vec(A) = vec(reshuffle(A)).
1946 V, I, J = [1]*mn, [], range(mn)
1947 for i in range(mn):
1948 (k, j) = divmod(i, m) # (j, k) are column-major matrix indices.
1949 indices = reindex((j, k), (m, n), dimensions)
1950 newIndices = tuple(indices[P[d]] for d in range(numDims))
1951 newDimensions = tuple(dimensions[P[d]] for d in range(numDims))
1952 (j, k) = reindex(newIndices, newDimensions, (m, n))
1953 I.append(k*m + j)
1954 R = cvxopt.spmatrix(V, I, J, (mn, mn), self._typecode)
1956 # Create the string.
1957 strArgs = [self.string, str(permutation).replace(" ", ""), order]
1959 if not guessDimensions:
1960 strArgs.insert(2, str(dimensions).replace(" ", ""))
1962 string = glyphs.shuffled(",".join(strArgs))
1964 # Finalize the new expression.
1965 shape = (m, n)
1966 coefs = {mtbs: R * coef for mtbs, coef in self._coefs.items()}
1968 return self._basetype(string, shape, coefs)
1970 @cached_property
1971 def sum(self):
1972 """Sum over all scalar elements of the expression."""
1973 # NOTE: glyphs.clever_dotp detects this case and uses the sum glyph.
1974 # NOTE: 1 on the left side in case self is complex.
1975 return (1 | self)
1977 @cached_property
1978 def rowsum(self):
1979 """Sum over the rows of the expression as a row vector."""
1980 from .algebra import J
1982 return J(self.shape[0]).T * self
1984 @cached_property
1985 def colsum(self):
1986 """Sum over the columns of the expression as a column vector."""
1987 from .algebra import J
1989 return self * J(self.shape[1])
1991 @cached_property
1992 def tr(self):
1993 """Trace of a square expression."""
1994 if not self.square:
1995 raise TypeError("Cannot compute the trace of non-square {}."
1996 .format(self.string))
1998 # NOTE: glyphs.clever_dotp detects this case and uses the trace glyph.
1999 # NOTE: "I" on the left side in case self is complex.
2000 return ("I" | self)
2002 def partial_trace(self, subsystems, dimensions=2):
2003 r"""Return the partial trace over selected subsystems.
2005 If the expression can be written as
2006 :math:`A_0 \otimes \cdots \otimes A_{n-1}` for matrices
2007 :math:`A_0, \ldots, A_{n-1}` with shapes given in ``dimensions``, then
2008 this returns :math:`B_0 \otimes \cdots \otimes B_{n-1}` with
2009 :math:`B_i = \operatorname{tr}(A_i)`, if ``i in subsystems`` (with
2010 :math:`i = -1` read as :math:`n-1`), and :math:`B_i = A_i`, otherwise.
2012 :param subsystems: A collection of or a single subystem number, indexed
2013 from zero, corresponding to subsystems that shall be traced over.
2014 The value :math:`-1` refers to the last subsystem.
2015 :type subsystems: int or tuple or list
2017 :param dimensions: Either an integer :math:`d` so that the subsystems
2018 are assumed to be all of shape :math:`d \times d`, or a sequence of
2019 subsystem shapes where an integer :math:`d` within the sequence is
2020 read as :math:`d \times d`. In any case, the elementwise product
2021 over all subsystem shapes must equal the expression's shape.
2022 :type dimensions: int or tuple or list
2024 :raises TypeError: If the subsystems do not match the expression or if
2025 a non-square subsystem is to be traced over.
2026 :raises IndexError: If the subsystem selection is invalid in any other
2027 way.
2029 :Example:
2031 >>> from picos import Constant
2032 >>> A = Constant("A", range(16), (4, 4))
2033 >>> print(A) #doctest: +NORMALIZE_WHITESPACE
2034 [ 0.00e+00 4.00e+00 8.00e+00 1.20e+01]
2035 [ 1.00e+00 5.00e+00 9.00e+00 1.30e+01]
2036 [ 2.00e+00 6.00e+00 1.00e+01 1.40e+01]
2037 [ 3.00e+00 7.00e+00 1.10e+01 1.50e+01]
2038 >>> A0 = A.partial_trace(0); A0
2039 <2×2 Real Constant: A.{tr([2×2])⊗[2×2]}>
2040 >>> print(A0) #doctest: +NORMALIZE_WHITESPACE
2041 [ 1.00e+01 1.80e+01]
2042 [ 1.20e+01 2.00e+01]
2043 >>> A1 = A.partial_trace(1); A1
2044 <2×2 Real Constant: A.{[2×2]⊗tr([2×2])}>
2045 >>> print(A1) #doctest: +NORMALIZE_WHITESPACE
2046 [ 5.00e+00 2.10e+01]
2047 [ 9.00e+00 2.50e+01]
2048 """
2049 # Shape of the original matrix.
2050 m, n = self._shape
2052 if isinstance(dimensions, int):
2053 dimensions = self._square_equal_subsystem_dims(dimensions)
2054 else:
2055 dimensions = [
2056 (d, d) if isinstance(d, int) else d for d in dimensions]
2058 if reduce(
2059 lambda x, y: (x[0]*y[0], x[1]*y[1]), dimensions) != self._shape:
2060 raise TypeError("Subsystem dimensions do not match expression.")
2062 if isinstance(subsystems, int):
2063 subsystems = (subsystems,)
2064 elif not subsystems:
2065 return self
2067 numSys = len(dimensions)
2068 subsystems = set(numSys - 1 if sys == -1 else sys for sys in subsystems)
2070 for sys in subsystems:
2071 if not isinstance(sys, int):
2072 raise IndexError("Subsystem indices must be integer, not {}."
2073 .format(type(sys).__name__))
2074 elif sys < 0:
2075 raise IndexError("Subsystem indices must be nonnegative.")
2076 elif sys >= numSys:
2077 raise IndexError(
2078 "Subsystem index {} out of range for {} systems total."
2079 .format(sys, numSys))
2080 elif dimensions[sys][0] != dimensions[sys][1]:
2081 raise TypeError(
2082 "Subsystem index {} refers to a non-square subsystem that "
2083 "cannot be traced over.".format(sys))
2085 # If all subsystems are traced over, this is the regular trace.
2086 if len(subsystems) == numSys:
2087 return self.tr
2089 # Prepare sparse T such that T·vec(A) = vec(partial_trace(A)).
2090 T = []
2092 # Compute factors of T, one for each subsystem being traced over.
2093 ibh, ibw = m, n
2094 sysStrings = None
2095 for sys in range(numSys):
2096 # Shape of current system.
2097 p, q = dimensions[sys]
2098 sysString = glyphs.matrix(glyphs.shape((p, q)))
2100 # Heigh/width of "outer" blocks whose relative position is
2101 # maintained but that are reduced independently to the size of an
2102 # inner block if the current system is to be traced over.
2103 obh, obw = ibh, ibw
2105 # Height/width of "inner" blocks that are summed if they are
2106 # main-diagonal blocks of an outer block.
2107 ibh, ibw = obh // p, obw // q
2109 # Only trace over selected subsystems.
2110 if sys not in subsystems:
2111 sysStrings = glyphs.kron(sysStrings, sysString) \
2112 if sysStrings else sysString
2113 continue
2114 else:
2115 sysStrings = glyphs.kron(sysStrings, glyphs.trace(sysString)) \
2116 if sysStrings else glyphs.trace(sysString)
2118 # Shape of new matrix.
2119 assert p == q
2120 mN, nN = m // p, n // p
2122 # Prepare one factor of T.
2123 oldLen = m * n
2124 newLen = mN * nN
2125 V, I, J = [1]*(newLen*p), [], []
2126 shape = (newLen, oldLen)
2128 # Determine the summands that make up each entry of the new matrix.
2129 for viN in range(newLen):
2130 # A column/row pair that represents a scalar in the new matrix
2131 # and a number of p scalars within different on-diagonal inner
2132 # blocks but within the same outer block in the old matrix.
2133 cN, rN = divmod(viN, mN)
2135 # Index pair (obi, obj) for outer block in question, row/column
2136 # pair (ibr, ibc) identifies the scalar within each inner block.
2137 obi, ibr = divmod(rN, ibh)
2138 obj, ibc = divmod(cN, ibw)
2140 # Collect summands for the current entry of the new matrix; one
2141 # scalar per on-diagonal inner block.
2142 for k in range(p):
2143 rO = obi*obh + k*ibh + ibr
2144 cO = obj*obw + k*ibw + ibc
2145 I.append(viN)
2146 J.append(cO*m + rO)
2148 # Store the operator that performs the current subsystem trace.
2149 T.insert(0, cvxopt.spmatrix(V, I, J, shape, self._typecode))
2151 # Make next iteration work on the new matrix.
2152 m, n = mN, nN
2154 # Multiply out the linear partial trace operator T.
2155 # TODO: Fast matrix multiplication dynamic program useful here?
2156 T = reduce(lambda A, B: A*B, T)
2158 string = glyphs.ptrace_(self.string, sysStrings)
2159 shape = (m, n)
2160 coefs = {mtbs: T * coef for mtbs, coef in self._coefs.items()}
2162 return self._basetype(string, shape, coefs)
2164 @cached_property
2165 def tr0(self):
2166 r"""Expression with the first :math:`2 \times 2` subsystem traced out.
2168 Only available for a :math:`2^k \times 2^k` matrix with all subsystems
2169 of shape :math:`2 \times 2`. Use :meth:`partial_trace` otherwise.
2170 """
2171 return self.partial_trace(subsystems=0)
2173 @cached_property
2174 def tr1(self):
2175 r"""Expression with the second :math:`2 \times 2` subsystem traced out.
2177 Only available for a :math:`2^k \times 2^k` matrix with all subsystems
2178 of shape :math:`2 \times 2`. Use :meth:`partial_trace` otherwise.
2179 """
2180 return self.partial_trace(subsystems=1)
2182 @cached_property
2183 def tr2(self):
2184 r"""Expression with the third :math:`2 \times 2` subsystem traced out.
2186 Only available for a :math:`2^k \times 2^k` matrix with all subsystems
2187 of shape :math:`2 \times 2`. Use :meth:`partial_trace` otherwise.
2188 """
2189 return self.partial_trace(subsystems=2)
2191 @cached_property
2192 def tr3(self):
2193 r"""Expression with the fourth :math:`2 \times 2` subsystem traced out.
2195 Only available for a :math:`2^k \times 2^k` matrix with all subsystems
2196 of shape :math:`2 \times 2`. Use :meth:`partial_trace` otherwise.
2197 """
2198 return self.partial_trace(subsystems=3)
2200 @cached_property
2201 def trl(self):
2202 r"""Expression with the last :math:`2 \times 2` subsystem traced out.
2204 Only available for a :math:`2^k \times 2^k` matrix with all subsystems
2205 of shape :math:`2 \times 2`. Use :meth:`partial_trace` otherwise.
2206 """
2207 return self.partial_trace(subsystems=-1)
2209 @cached_property
2210 def vec(self):
2211 """Column-major vectorization of the expression as a column vector.
2213 .. note::
2214 Given an expression ``A``, ``A.vec`` and ``A[:]`` produce the same
2215 result (up to its string description) but ``A.vec`` is faster and
2216 its result is cached.
2218 :Example:
2220 >>> from picos import Constant
2221 >>> A = Constant("A", [[1, 2], [3, 4]])
2222 >>> A.vec.equals(A[:])
2223 True
2224 >>> A[:] is A[:]
2225 False
2226 >>> A.vec is A.vec
2227 True
2228 """
2229 if self._shape[1] == 1:
2230 return self
2231 else:
2232 return self._basetype(
2233 glyphs.vec(self.string), (len(self), 1), self._coefs)
2235 def dupvec(self, n):
2236 """Return a (repeated) column-major vectorization of the expression.
2238 :param int n: Number of times to duplicate the vectorization.
2240 :returns: A column vector.
2242 :Example:
2244 >>> from picos import Constant
2245 >>> A = Constant("A", [[1, 2], [3, 4]])
2246 >>> A.dupvec(1) is A.vec
2247 True
2248 >>> A.dupvec(3).equals(A.vec // A.vec // A.vec)
2249 True
2250 """
2251 if not isinstance(n, int):
2252 raise TypeError("Number of copies must be integer.")
2254 if n < 1:
2255 raise ValueError("Number of copies must be positive.")
2257 if n == 1:
2258 return self.vec
2259 else:
2260 string = glyphs.vec(glyphs.comma(self.string, n))
2261 shape = (len(self)*n, 1)
2262 coefs = {mtbs: cvxopt_vcat([coef]*n)
2263 for mtbs, coef in self._coefs.items()}
2265 return self._basetype(string, shape, coefs)
2267 @cached_property
2268 def trilvec(self):
2269 r"""Column-major vectorization of the lower triangular part.
2271 :returns:
2272 A column vector of all elements :math:`A_{ij}` that satisfy
2273 :math:`i \geq j`.
2275 .. note::
2277 If you want a row-major vectorization instead, write ``A.T.triuvec``
2278 instead of ``A.trilvec``.
2280 :Example:
2282 >>> from picos import Constant
2283 >>> A = Constant("A", [[1, 2], [3, 4], [5, 6]])
2284 >>> print(A)
2285 [ 1.00e+00 2.00e+00]
2286 [ 3.00e+00 4.00e+00]
2287 [ 5.00e+00 6.00e+00]
2288 >>> print(A.trilvec)
2289 [ 1.00e+00]
2290 [ 3.00e+00]
2291 [ 5.00e+00]
2292 [ 4.00e+00]
2293 [ 6.00e+00]
2294 """
2295 m, n = self._shape
2297 if n == 1: # Column vector or scalar.
2298 return self
2299 elif m == 1: # Row vector.
2300 return self[0]
2302 # Build a transformation D such that D·vec(A) = trilvec(A).
2303 rows = [j*m + i for j in range(n) for i in range(m) if i >= j]
2304 d = len(rows)
2305 D = cvxopt.spmatrix([1]*d, range(d), rows, (d, len(self)))
2307 string = glyphs.trilvec(self.string)
2308 shape = (d, 1)
2309 coefs = {mtbs: D*coef for mtbs, coef in self._coefs.items()}
2311 return self._basetype(string, shape, coefs)
2313 @cached_property
2314 def triuvec(self):
2315 r"""Column-major vectorization of the upper triangular part.
2317 :returns:
2318 A column vector of all elements :math:`A_{ij}` that satisfy
2319 :math:`i \leq j`.
2321 .. note::
2323 If you want a row-major vectorization instead, write ``A.T.trilvec``
2324 instead of ``A.triuvec``.
2326 :Example:
2328 >>> from picos import Constant
2329 >>> A = Constant("A", [[1, 2, 3], [4, 5, 6]])
2330 >>> print(A)
2331 [ 1.00e+00 2.00e+00 3.00e+00]
2332 [ 4.00e+00 5.00e+00 6.00e+00]
2333 >>> print(A.triuvec)
2334 [ 1.00e+00]
2335 [ 2.00e+00]
2336 [ 5.00e+00]
2337 [ 3.00e+00]
2338 [ 6.00e+00]
2339 """
2340 m, n = self._shape
2342 if m == 1: # Row vector or scalar.
2343 return self
2344 elif n == 1: # Column vector.
2345 return self[0]
2347 # Build a transformation D such that D·vec(A) = triuvec(A).
2348 rows = [j*m + i for j in range(n) for i in range(m) if i <= j]
2349 d = len(rows)
2350 D = cvxopt.spmatrix([1]*d, range(d), rows, (d, len(self)))
2352 string = glyphs.triuvec(self.string)
2353 shape = (d, 1)
2354 coefs = {mtbs: D*coef for mtbs, coef in self._coefs.items()}
2356 return self._basetype(string, shape, coefs)
2358 @cached_property
2359 def svec(self):
2360 """An isometric vectorization of a symmetric or Hermitian expression.
2362 In the real symmetric case
2364 - the vectorization format is precisely the one define in [svec]_,
2365 - the vectorization is isometric and isomorphic, and
2366 - this is the same vectorization as used internally by the
2367 :class:`~.variables.SymmetricVariable` class.
2369 In the complex hermitian case
2371 - the same format is used, now resulting in a complex vector,
2372 - the vectorization is isometric but **not** isomorphic as there are
2373 guaranteed zeros in the imaginary part of the vector, and
2374 - this is **not** the same vectorization as the isomorphic,
2375 real-valued one used by :class:`~.variables.HermitianVariable`.
2377 The reverse operation is denoted by :attr:`desvec` in either case.
2379 :raises TypeError:
2380 If the expression is not square.
2382 :raises ValueError:
2383 If the expression is not hermitian.
2384 """
2385 if not self.square:
2386 raise TypeError("Can only compute svec for a square matrix, not for"
2387 " the {} expression {}.".format(self._shape, self.string))
2388 elif not self.hermitian:
2389 raise ValueError("Cannot compute svec for the non-hermitian "
2390 "expression {}.".format(self.string))
2392 vec = SymmetricVectorization(self._shape)
2393 V = vec._full2special
2395 string = glyphs.svec(self.string)
2397 if self.isreal:
2398 vec = SymmetricVectorization(self._shape)
2399 V = vec._full2special
2400 coefs = {mtbs: V*coef for mtbs, coef in self._coefs.items()}
2401 result = self._basetype(string, (vec.dim, 1), coefs)
2402 else:
2403 # NOTE: We need to reproduce svec using triuvec because, for numeric
2404 # reasons, SymmetricVectorization averages the elements in the
2405 # lower and upper triangular part. For symmetric matrices,
2406 # this is equivalent to the formal definition of svec found in
2407 # Datorro's book. For hermitian matrices however it is not.
2408 real_svec = self.real.svec
2409 imag_svec = 2**0.5 * 1j * self.imag.triuvec
2411 result = (real_svec + imag_svec).renamed(string)
2413 with unlocked_cached_properties(result):
2414 result.desvec = self
2416 return result
2418 @cached_property
2419 def desvec(self):
2420 r"""The reverse operation of :attr:`svec`.
2422 :raises TypeError:
2423 If the expression is not a vector or has a dimension outside the
2424 permissible set
2426 .. math::
2428 \left\{ \frac{n(n + 1)}{2}
2429 \mid n \in \mathbb{Z}_{\geq 1} \right\}
2430 = \left\{ n \in \mathbb{Z}_{\geq 1}
2431 \mid \frac{1}{2} \left( \sqrt{8n + 1} - 1 \right)
2432 \in \mathbb{Z}_{\geq 1} \right\}.
2434 :raises ValueError:
2435 In the case of a complex vector, If the vector is not in the image
2436 of :attr:`svec`.
2437 """
2438 if 1 not in self.shape:
2439 raise TypeError(
2440 "Can only compute desvec for a vector, not for the {} "
2441 "expression {}.".format(glyphs.shape(self._shape), self.string))
2443 n = 0.5*((8*len(self) + 1)**0.5 - 1)
2444 if int(n) != n:
2445 raise TypeError("Cannot compute desvec for the {}-dimensional "
2446 "vector {} as this size is not a possible outcome of svec."
2447 .format(len(self), self.string))
2448 n = int(n)
2450 vec = SymmetricVectorization((n, n))
2451 D = vec._special2full
2452 string = glyphs.desvec(self.string)
2454 if self.isreal:
2455 coefs = {mtbs: D*coef for mtbs, coef in self._coefs.items()}
2456 result = self._basetype(string, (n, n), coefs)
2458 assert result.hermitian
2459 else:
2460 # NOTE: While :attr:`svec` performs essentially the same operation
2461 # for both symmetric and hermitian matrices, we now need to
2462 # treat the imaginary separately as the imaginary part of a
2463 # hermitian matrix is skew-symmetric instead of symmetric.
2464 V, I, J = [], D.I, D.J
2465 for v, i, j in zip(D.V, I, J):
2466 V.append(1j*v if i % n <= i // n else -1j*v)
2467 C = cvxopt.spmatrix(V, I, J, D.size, tc="z")
2469 real_desvec = self.real.desvec
2470 imag_desvec = self._basetype(string, (n, n), {
2471 mtbs: C*coef for mtbs, coef in self.imag._coefs.items()})
2473 result = (real_desvec + imag_desvec).renamed(string)
2475 if not result.hermitian:
2476 raise ValueError("The complex vector {} is not in the image of "
2477 "svec. Note that svec is not bijective in the complex case."
2478 .format(self.string))
2480 with unlocked_cached_properties(result):
2481 result.svec = self
2483 return result
2485 def dupdiag(self, n):
2486 """Return a matrix with the (repeated) expression on the diagonal.
2488 Vectorization is performed in column-major order.
2490 :param int n: Number of times to duplicate the vectorization.
2491 """
2492 from .algebra import I
2494 if self.scalar:
2495 return self * I(n)
2497 # Vectorize and duplicate the expression.
2498 vec = self.dupvec(n)
2499 d = len(vec)
2501 # Build a transformation D such that D·vec(A) = vec(diag(vec(A))).
2502 ones = [1]*d
2503 D = cvxopt.spdiag(ones)[:]
2504 D = cvxopt.spmatrix(ones, D.I, range(d), (D.size[0], d))
2506 string = glyphs.diag(vec.string)
2507 shape = (d, d)
2508 coefs = {mtbs: D*coef for mtbs, coef in vec._coefs.items()}
2510 return self._basetype(string, shape, coefs)
2512 @cached_property
2513 def diag(self):
2514 """Diagonal matrix with the expression on the main diagonal.
2516 Vectorization is performed in column-major order.
2517 """
2518 from .algebra import O, I
2520 if self.is0:
2521 return O(len(self), len(self))
2522 elif self.is1:
2523 return I(len(self))
2524 else:
2525 return self.dupdiag(1)
2527 @cached_property
2528 def maindiag(self):
2529 """The main diagonal of the expression as a column vector."""
2530 if 1 in self._shape:
2531 return self[0]
2533 # Build a transformation D such that D·vec(A) = diag(A).
2534 step = self._shape[0] + 1
2535 rows = [i*step for i in range(min(self._shape))]
2536 d = len(rows)
2537 D = cvxopt.spmatrix([1]*d, range(d), rows, (d, len(self)))
2539 string = glyphs.maindiag(self.string)
2540 shape = (d, 1)
2541 coefs = {mtbs: D*coef for mtbs, coef in self._coefs.items()}
2543 return self._basetype(string, shape, coefs)
2545 def factor_out(self, mutable):
2546 r"""Factor out a single mutable from a vector expression.
2548 If this expression is a column vector :math:`a` that depends on some
2549 mutable :math:`x` with a trivial internal vectorization format (i.e.
2550 :class:`~.vectorizations.FullVectorization`), then this method, called
2551 with :math:`x` as its argument, returns a pair of expressions
2552 :math:`(a_x, a_0)` such that :math:`a = a_x\operatorname{vec}(x) + a_0`.
2554 :returns:
2555 Two refined :class:`BiaffineExpression` instances that do not depend
2556 on ``mutable``.
2558 :raises TypeError:
2559 If the expression is not a column vector or if ``mutable`` is not a
2560 :class:`~.mutable.Mutable` or does not have a trivial vectorization
2561 format.
2563 :raises LookupError:
2564 If the expression does not depend on ``mutable``.
2566 :Example:
2568 >>> from picos import RealVariable
2569 >>> from picos.uncertain import UnitBallPerturbationSet
2570 >>> x = RealVariable("x", 3)
2571 >>> z = UnitBallPerturbationSet("z", 3).parameter
2572 >>> a = ((2*x + 3)^z) + 4*x + 5; a
2573 <3×1 Uncertain Affine Expression: (2·x + [3])⊙z + 4·x + [5]>
2574 >>> sorted(a.mutables, key=lambda mtb: mtb.name)
2575 [<3×1 Real Variable: x>, <3×1 Perturbation: z>]
2576 >>> az, a0 = a.factor_out(z)
2577 >>> az
2578 <3×3 Real Affine Expression: ((2·x + [3])⊙z + 4·x + [5])_z>
2579 >>> a0
2580 <3×1 Real Affine Expression: ((2·x + [3])⊙z + 4·x + [5])_0>
2581 >>> sorted(az.mutables.union(a0.mutables), key=lambda mtb: mtb.name)
2582 [<3×1 Real Variable: x>]
2583 >>> (az*z + a0).equals(a)
2584 True
2585 """
2586 if self._shape[1] != 1:
2587 raise TypeError(
2588 "Can only factor out mutables from column vectors but {} has "
2589 "a shape of {}.".format(self.string, glyphs.shape(self._shape)))
2591 mtb = mutable
2593 if not isinstance(mtb, Mutable):
2594 raise TypeError("Can only factor out mutables, not instances of {}."
2595 .format(type(mtb).__name__))
2597 if not isinstance(mtb._vec, FullVectorization):
2598 raise TypeError("Can only factor out mutables with a trivial "
2599 "vectorization format but {} uses {}."
2600 .format(mtb.name, type(mtb._vec).__name__))
2602 if mtb not in self.mutables:
2603 raise LookupError("Cannot factor out {} from {} as the latter does "
2604 "not depend on the former.".format(mtb.name, self.string))
2606 ax_string = glyphs.index(self.string, mtb.name)
2607 ax_shape = (self._shape[0], mtb.dim)
2609 # Recipe in "Robust conic optimization in Python" (Stahlberg 2020).
2610 ax_coefs = {}
2611 for mtbs, coef in self._coefs.items():
2612 if mtb not in set(mtbs):
2613 continue
2614 elif len(mtbs) == 1:
2615 assert mtbs[0] is mtb
2616 self._update_coefs(ax_coefs, (), coef[:])
2617 else:
2618 if mtbs[0] is mtb:
2619 other_mtb = mtbs[1]
2620 C = coef*cvxopt_K(other_mtb.dim, mtb.dim)
2621 else:
2622 assert mtbs[1] is mtb
2623 other_mtb = mtbs[0]
2624 C = coef
2626 d = other_mtb.dim
2627 C = cvxopt.sparse([C[:, i*d:(i+1)*d] for i in range(mtb.dim)])
2628 self._update_coefs(ax_coefs, (other_mtb,), C)
2630 ax = self._basetype(ax_string, ax_shape, ax_coefs)
2632 a0_string = glyphs.index(self.string, 0)
2633 a0_shape = self._shape
2634 a0_coefs = {M: C for M, C in self._coefs.items() if mtb not in set(M)}
2635 a0 = self._basetype(a0_string, a0_shape, a0_coefs)
2637 return ax.refined, a0.refined
2639 # --------------------------------------------------------------------------
2640 # Backwards compatibility methods.
2641 # --------------------------------------------------------------------------
2643 @classmethod
2644 @deprecated("2.0", useInstead="from_constant")
2645 def fromScalar(cls, scalar):
2646 """Create a class instance from a numeric scalar."""
2647 return cls.from_constant(scalar, (1, 1))
2649 @classmethod
2650 @deprecated("2.0", useInstead="from_constant")
2651 def fromMatrix(cls, matrix, size=None):
2652 """Create a class instance from a numeric matrix."""
2653 return cls.from_constant(matrix, size)
2655 @deprecated("2.0", useInstead="object.__xor__")
2656 def hadamard(self, fact):
2657 """Denote the elementwise (or Hadamard) product."""
2658 return self.__xor__(fact)
2660 @deprecated("2.0", useInstead="~.expression.Expression.constant")
2661 def isconstant(self):
2662 """Whether the expression involves no mutables."""
2663 return self.constant
2665 @deprecated("2.0", useInstead="equals")
2666 def same_as(self, other):
2667 """Check mathematical equality with another affine expression."""
2668 return self.equals(other)
2670 @deprecated("2.0", useInstead="T")
2671 def transpose(self):
2672 """Return the matrix transpose."""
2673 return self.T
2675 @cached_property
2676 @deprecated("2.0", useInstead="partial_transpose", decoratorLevel=1)
2677 def Tx(self):
2678 """Auto-detect few subsystems of same shape and transpose the last."""
2679 m, n = self._shape
2680 dims = None
2682 for k in range(2, int(math.log(min(m, n), 2)) + 1):
2683 p, q = int(round(m**(1.0/k))), int(round(n**(1.0/k)))
2684 if m == p**k and n == q**k:
2685 dims = ((p, q),)*k
2686 break
2688 if dims:
2689 return self.partial_transpose(subsystems=-1, dimensions=dims)
2690 else:
2691 raise RuntimeError("Failed to auto-detect subsystem dimensions for "
2692 "partial transposition: Only supported for {} matrices, {}."
2693 .format(glyphs.shape(
2694 (glyphs.power("m", "k"), glyphs.power("n", "k"))),
2695 glyphs.ge("k", 2)))
2697 @deprecated("2.0", useInstead="conj")
2698 def conjugate(self):
2699 """Return the complex conjugate."""
2700 return self.conj
2702 @deprecated("2.0", useInstead="H")
2703 def Htranspose(self):
2704 """Return the conjugate (or Hermitian) transpose."""
2705 return self.H
2707 @deprecated("2.0", reason="PICOS expressions are now immutable.")
2708 def copy(self):
2709 """Return a deep copy of the expression."""
2710 from copy import copy as cp
2712 return self._basetype(cp(self._symbStr), self._shape,
2713 {mtbs: cp(coef) for mtbs, coef in self._coefs.items()})
2715 @deprecated("2.0", reason="PICOS expressions are now immutable.")
2716 def soft_copy(self):
2717 """Return a shallow copy of the expression."""
2718 return self._basetype(self._symbStr, self._shape, self._coefs)
2721# --------------------------------------
2722__all__ = api_end(_API_START, globals())