Coverage for picos/expressions/exp_biaffine.py : 82.12%
Hot-keys on this page
r m x p toggle line displays
j k next/prev highlighted chunk
0 (zero) top of page
1 (one) first highlighted chunk
1# ------------------------------------------------------------------------------
2# Copyright (C) 2019-2020 Maximilian Stahlberg
3# Based on the original picos.expressions module by Guillaume Sagnol.
4#
5# This file is part of PICOS.
6#
7# PICOS is free software: you can redistribute it and/or modify it under the
8# terms of the GNU General Public License as published by the Free Software
9# Foundation, either version 3 of the License, or (at your option) any later
10# version.
11#
12# PICOS is distributed in the hope that it will be useful, but WITHOUT ANY
13# WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR
14# A PARTICULAR PURPOSE. See the GNU General Public License for more details.
15#
16# You should have received a copy of the GNU General Public License along with
17# this program. If not, see <http://www.gnu.org/licenses/>.
18# ------------------------------------------------------------------------------
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) if shape == (1, 1) else glyphs.matrix(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 (re-)cast to sparse matrices."""
367 return {mtbs: cvxopt.sparse(coef) for mtbs, coef in self._coefs.items()}
369 @cached_property
370 def _sparse_linear_coefs(self):
371 """Linear part coefficients cast to sparse and indexed by mutables."""
372 return {v[0]: c for v, c in self._sparse_coefs.items() if len(v) == 1}
374 def _reglyphed_string(self, mutable_name_map):
375 """The symbolic string with mutable names replaced."""
376 string = self.string
378 if isinstance(string, glyphs.GlStr):
379 string = string.reglyphed(mutable_name_map)
380 elif string in mutable_name_map:
381 # Handle corner cases like x + 0 for a mutable x which is not a
382 # mutable any more but still has a literal string name.
383 string = mutable_name_map[string]
385 return string
387 # --------------------------------------------------------------------------
388 # Abstract method implementations and overridings for Expression.
389 # --------------------------------------------------------------------------
391 def _get_value(self):
392 # Create a copy of the constant term.
393 value = self._constant_coef[:]
395 # Add value of the linear part.
396 for mtb, coef in self._linear_coefs.items():
397 summand = coef * mtb._get_internal_value()
399 if type(value) == type(summand):
400 value += summand
401 else:
402 # Exactly one of the matrices is sparse.
403 value = value + summand
405 # Add value of the bilinear part.
406 for (x, y), coef in self._bilinear_coefs.items():
407 xValue = x._get_internal_value()
408 yValue = y._get_internal_value()
409 summand = coef * (xValue*yValue.T)[:]
411 if type(value) == type(summand):
412 value += summand
413 else:
414 # Exactly one of the matrices is sparse.
415 value = value + summand
417 # Resize the value to the proper shape.
418 value.size = self.shape
420 return value
422 def _get_shape(self):
423 return self._shape
425 @cached_unary_operator
426 def _get_mutables(self):
427 return frozenset(mtb for mtbs in self._coefs for mtb in mtbs)
429 def _is_convex(self):
430 return not self._bilinear_coefs
432 def _is_concave(self):
433 return not self._bilinear_coefs
435 def _replace_mutables(self, mapping):
436 # Handle the base case where the affine expression is a mutable.
437 if self in mapping:
438 return mapping[self]
440 name_map = {old.name: new.string for old, new in mapping.items()}
441 string = self._reglyphed_string(name_map)
443 if all(isinstance(new, Mutable) for new in mapping.values()):
444 # Fast implementation for the basic case.
445 coefs = {tuple(mapping[mtb] for mtb in mtbs): coef
446 for mtbs, coef in self._coefs.items()}
447 else:
448 # Turn full mapping into an effective mapping.
449 mapping = {o: n for o, n in mapping.items() if o is not n}
451 coefs = {}
452 for old_mtbs, old_coef in self._coefs.items():
453 if not any(old_mtb in mapping for old_mtb in old_mtbs):
454 self._update_coefs(coefs, old_mtbs, old_coef)
455 elif len(old_mtbs) == 1:
456 assert old_mtbs[0] in mapping
457 new = mapping[old_mtbs[0]]
458 for new_mtb, new_coef in new._linear_coefs.items():
459 self._update_coefs(coefs, (new_mtb,), old_coef*new_coef)
460 self._update_coefs(coefs, (), old_coef*new._constant_coef)
461 elif old_mtbs[0] in mapping and old_mtbs[1] not in mapping:
462 new = mapping[old_mtbs[0]]
463 old_mtb = old_mtbs[1]
464 for new_mtb, new_coef in new._linear_coefs.items():
465 self._update_coefs(coefs, (new_mtb, old_mtb),
466 old_coef*(left_kronecker_I(new_coef, old_mtb.dim)))
467 self._update_coefs(coefs, (old_mtb,), old_coef*(
468 left_kronecker_I(new._constant_coef, old_mtb.dim)))
469 elif old_mtbs[0] not in mapping and old_mtbs[1] in mapping:
470 new = mapping[old_mtbs[1]]
471 old_mtb = old_mtbs[0]
472 for new_mtb, new_coef in new._linear_coefs.items():
473 self._update_coefs(coefs, (old_mtb, new_mtb),
474 old_coef*(right_kronecker_I(new_coef, old_mtb.dim)))
475 self._update_coefs(coefs, (old_mtb,), old_coef*(
476 right_kronecker_I(new._constant_coef, old_mtb.dim)))
477 elif old_mtbs[0] in mapping and old_mtbs[1] in mapping:
478 new1 = mapping[old_mtbs[0]]
479 new2 = mapping[old_mtbs[1]]
480 if isinstance(new1, Mutable) and isinstance(new2, Mutable):
481 self._update_coefs(coefs, (new1, new2), old_coef)
482 else:
483 raise NotImplementedError(
484 "Replacing both mutables in a bilinear term is not "
485 "supported unless both are replaced with mutables. "
486 "The effective mapping was: {}".format(mapping))
487 else:
488 assert False
490 return self._basetype(string, self._shape, coefs)
492 def _freeze_mutables(self, freeze):
493 string = self._reglyphed_string(
494 {mtb.name: glyphs.frozen(mtb.name) for mtb in freeze})
496 coefs = {}
497 for mtbs, coef in self._coefs.items():
498 if not any(mtb in freeze for mtb in mtbs):
499 self._update_coefs(coefs, mtbs, coef)
500 elif len(mtbs) == 1:
501 assert mtbs[0] in freeze
502 self._update_coefs(coefs, (), coef*mtbs[0].internal_value)
503 elif mtbs[0] in freeze and mtbs[1] in freeze:
504 C = coef*(mtbs[0].internal_value*mtbs[1].internal_value.T)[:]
505 self._update_coefs(coefs, (), C)
506 elif mtbs[0] in freeze and mtbs[1] not in freeze:
507 C = coef*left_kronecker_I(mtbs[0].internal_value, mtbs[1].dim)
508 self._update_coefs(coefs, (mtbs[1],), C)
509 elif mtbs[0] not in freeze and mtbs[1] in freeze:
510 C = coef*right_kronecker_I(mtbs[1].internal_value, mtbs[0].dim)
511 self._update_coefs(coefs, (mtbs[0],), C)
512 else:
513 assert False
515 return self._basetype(string, self._shape, coefs)
517 # --------------------------------------------------------------------------
518 # Python special method implementations.
519 # --------------------------------------------------------------------------
521 def __len__(self):
522 # Faster version that overrides Expression.__len__.
523 return self._shape[0] * self._shape[1]
525 def __getitem__(self, index):
526 def slice2range(s, length):
527 """Transform a :class:`slice` to a :class:`range`."""
528 assert isinstance(s, slice)
530 # Plug in slice's default values.
531 ss = s.step if s.step else 1
532 if ss > 0:
533 sa = s.start if s.start is not None else 0
534 sb = s.stop if s.stop is not None else length
535 else:
536 assert ss < 0
537 sa = s.start if s.start is not None else length - 1
538 sb = s.stop # Keep it None as -1 would mean length - 1.
540 # Wrap negative indices (once).
541 ra = length + sa if sa < 0 else sa
542 if sb is None:
543 # This is the only case where we give a negative index to range.
544 rb = -1
545 else:
546 rb = length + sb if sb < 0 else sb
548 # Clamp out-of-bound indices.
549 ra = min(max(0, ra), length - 1)
550 rb = min(max(-1, rb), length)
552 r = range(ra, rb, ss)
554 if not r:
555 raise IndexError("Empty slice.")
557 return r
559 def range2slice(r, length):
560 """Transform a :class:`range` to a :class:`slice`, if possible.
562 :raises ValueError: If the input cannot be expressed as a slice.
563 """
564 assert isinstance(r, range)
566 if not r:
567 raise IndexError("Empty range.")
569 ra = r.start
570 rb = r.stop
571 rs = r.step
573 if rs > 0:
574 if ra < 0 or rb > length:
575 raise ValueError(
576 "Out-of-bounds range cannot be represented as a slice.")
577 else:
578 assert rs < 0
579 if ra >= length or rb < -1:
580 raise ValueError(
581 "Out-of-bounds range cannot be represented as a slice.")
583 if rb == -1:
584 rb = None
586 return slice(ra, rb, rs)
588 def list2slice(l, length):
589 """Transform a :class:`list` to a :class:`slice`, if possible.
591 :raises TypeError: If the input is not an integer sequence.
592 :raises ValueError: If the input cannot be expressed as a slice.
593 """
594 return range2slice(detect_range(l), length)
596 def slice2str(s):
597 """Return the short string that produced a :class:`slice`."""
598 assert isinstance(s, slice)
600 startStr = str(s.start) if s.start is not None else ""
601 stopStr = str(s.stop) if s.stop is not None else ""
602 if s.step in (None, 1):
603 if s.start is not None and s.stop is not None \
604 and s.stop == s.start + 1:
605 return startStr
606 else:
607 return "{}:{}".format(startStr, stopStr)
608 else:
609 return "{}:{}:{}".format(startStr, stopStr, str(s.step))
611 def list2str(l, length):
612 """Return a short string represnetation of a :class:`list`."""
613 assert isinstance(l, list)
615 # Extract integers wrapped in a list.
616 if len(l) == 1:
617 return str(l[0])
619 # Try to convert the list to a slice.
620 try:
621 l = list2slice(l, length)
622 except (ValueError, RuntimeError):
623 pass
625 if isinstance(l, list):
626 if len(l) > 4:
627 return glyphs.shortint(l[0], l[-1])
628 else:
629 return str(l).replace(" ", "")
630 else:
631 return slice2str(l)
633 def any2str(a, length):
634 if isinstance(a, slice):
635 return slice2str(a)
636 elif isinstance(a, list):
637 return list2str(a, length)
638 else:
639 assert False
641 m, n = self._shape
642 indexStr = None
643 isIntList = False
645 # Turn the index expression into a mutable list of one index per axis.
646 if isinstance(index, tuple): # Multiple axis slicing.
647 index = list(index)
648 elif isinstance(index, dict): # Arbitrary matrix element selection.
649 if len(index) != 2:
650 raise TypeError("When slicing with a dictionary, there must be "
651 "exactly two keys.")
653 try:
654 i, j = sorted(index.keys())
655 except TypeError as error:
656 raise TypeError("When slicing with a dictionary, the two keys "
657 "must be comparable.") from error
659 I, J = index[i], index[j]
661 try:
662 I = load_dense_data(I, typecode="i", alwaysCopy=False)[0]
663 J = load_dense_data(J, typecode="i", alwaysCopy=False)[0]
665 if 1 not in I.size or 1 not in J.size:
666 raise TypeError("At least one of the objects is not flat "
667 "but a proper matrix.")
669 if len(I) != len(J):
670 raise TypeError("The objects do not have the same length.")
672 I, J = list(I), list(J)
673 except Exception as error:
674 raise TypeError("Loading a sparse index vector pair for {} from"
675 " objects of type {} and {} failed: {}".format(self.string,
676 type(index[i]).__name__, type(index[j]).__name__, error)) \
677 from None
679 # Represent the selection as a global index list.
680 index = [[i + j*m for i, j in zip(I, J)]]
682 # Use a special index string.
683 indexStr = glyphs.size(list2str(I, n), list2str(J, m))
685 # Don't invoke load_dense_data on the list.
686 isIntList = True
687 else: # Global indexing.
688 index = [index]
690 # Make sure either global or row/column indexing is used.
691 if not index:
692 raise IndexError("Empty index.")
693 elif len(index) > 2:
694 raise IndexError(
695 "PICOS expressions do not have a third axis to slice.")
697 # Turn the indices for each axis into either a slice or a list.
698 for axis, selection in enumerate(index):
699 # Convert anything that is not a slice, including scalars and lists
700 # that are not confirmed integer, to an integer row or column
701 # vector, then (back) to a list.
702 if not isIntList and not isinstance(selection, slice):
703 try:
704 matrix = load_dense_data(
705 selection, typecode="i", alwaysCopy=False)[0]
707 if 1 not in matrix.size:
708 raise TypeError("The object is not flat but a {} shaped"
709 " matrix.".format(glyphs.shape(matrix.size)))
711 selection = list(matrix)
712 except Exception as error:
713 raise TypeError("Loading a slicing index vector for axis {}"
714 " of {} from an object of type {} failed: {}".format(
715 axis, self.string, type(selection).__name__, error)) \
716 from None
718 index[axis] = selection
720 # Build index string, retrieve new shape, finalize index.
721 if len(index) == 1: # Global indexing.
722 index = index[0]
724 if isinstance(index, slice):
725 shape = len(slice2range(index, len(self)))
726 else:
727 shape = len(index)
729 if indexStr is None:
730 indexStr = any2str(index, len(self))
731 else: # Multiple axis slicing.
732 if indexStr is None:
733 indexStr = "{},{}".format(
734 any2str(index[0], m), any2str(index[1], n))
736 # Convert to a global index list.
737 RC, shape = [], []
738 for axis, selection in enumerate(index):
739 k = self._shape[axis]
741 if isinstance(selection, slice):
742 # Turn the slice into an iterable range.
743 selection = slice2range(selection, k)
745 # All indices from a slice are nonnegative.
746 assert all(i >= 0 for i in selection)
748 if isinstance(selection, list):
749 # Wrap once. This is consistent with CVXOPT.
750 selection = [i if i >= 0 else k + i for i in selection]
752 # Perform a partial out-of-bounds check.
753 if any(i < 0 for i in selection):
754 raise IndexError(
755 "Out-of-bounds access along axis {}.".format(axis))
757 # Complete the check for out-of-bounds access.
758 if any(i >= k for i in selection):
759 raise IndexError(
760 "Out-of-bounds access along axis {}.".format(axis))
762 RC.append(selection)
763 shape.append(len(selection))
765 rows, cols = RC
766 index = [i + j*m for j in cols for i in rows]
768 # Finalize the string.
769 string = glyphs.slice(self.string, indexStr)
771 # Retrieve new coefficients and constant term.
772 coefs = {mtbs: coef[index, :] for mtbs, coef in self._coefs.items()}
774 return self._basetype(string, shape, coefs)
776 @convert_operands(sameShape=True)
777 @refine_operands(stop_at_affine=True)
778 def __add__(self, other):
779 """Denote addition from the right hand side."""
780 if not isinstance(other, BiaffineExpression):
781 return NotImplemented
783 string = glyphs.clever_add(self.string, other.string)
785 coefs = {}
786 for mtbs, coef in self._coefs.items():
787 coefs[mtbs] = coef + other._coefs[mtbs] \
788 if mtbs in other._coefs else coef
789 for mtbs, coef in other._coefs.items():
790 coefs.setdefault(mtbs, coef)
792 return self._common_basetype(other)(string, self._shape, coefs)
794 @convert_operands(sameShape=True)
795 @refine_operands(stop_at_affine=True)
796 def __radd__(self, other):
797 """Denote addition from the left hand side."""
798 if not isinstance(other, BiaffineExpression):
799 return NotImplemented
801 return other.__add__(self)
803 @convert_operands(sameShape=True)
804 @refine_operands(stop_at_affine=True)
805 def __sub__(self, other):
806 """Denote substraction from the right hand side."""
807 if not isinstance(other, BiaffineExpression):
808 return NotImplemented
810 string = glyphs.clever_sub(self.string, other.string)
812 coefs = {}
813 for mtbs, coef in self._coefs.items():
814 coefs[mtbs] = coef - other._coefs[mtbs] \
815 if mtbs in other._coefs else coef
816 for mtbs, coef in other._coefs.items():
817 coefs.setdefault(mtbs, -coef)
819 return self._common_basetype(other)(string, self._shape, coefs)
821 @convert_operands(sameShape=True)
822 @refine_operands(stop_at_affine=True)
823 def __rsub__(self, other):
824 """Denote substraction with self on the right hand side."""
825 if not isinstance(other, BiaffineExpression):
826 return NotImplemented
828 return other.__sub__(self)
830 def __pos__(self):
831 """Denote identity."""
832 return self
834 @cached_selfinverse_unary_operator
835 def __neg__(self):
836 """Denote negation."""
837 string = glyphs.clever_neg(self.string)
838 coefs = {mtbs: -coef for mtbs, coef in self._coefs.items()}
840 return self._basetype(string, self._shape, coefs)
842 @convert_operands(rMatMul=True)
843 @refine_operands(stop_at_affine=True)
844 def __mul__(self, other):
845 """Denote matrix multiplication from the right hand side."""
846 if not isinstance(other, BiaffineExpression):
847 return NotImplemented
849 string = glyphs.clever_mul(self.string, other.string)
851 m, n = self._shape
852 p, q = other._shape
854 # Handle or prepare multiplication with a scalar.
855 if 1 in (m*n, p*q):
856 if self.constant or other.constant:
857 # Handle all cases involving a constant operand immediately.
858 if self.constant:
859 lhs = self._constant_coef
860 RHS = other._coefs
861 else:
862 lhs = other._constant_coef
863 RHS = self._coefs
865 shape = other._shape if m*n == 1 else self._shape
867 # Work around CVXOPT not broadcasting sparse scalars.
868 if lhs.size == (1, 1):
869 lhs = lhs[0]
871 return self._common_basetype(other)(
872 string, shape, {mtbs: lhs*rhs for mtbs, rhs in RHS.items()})
873 elif n == p:
874 # Regular matrix multiplication already works.
875 pass
876 elif m*n == 1:
877 # Prepare a regular matrix multiplication.
878 # HACK: Replaces self.
879 self = self*cvxopt.spdiag([1]*p)
880 m, n = self._shape
881 else:
882 # Prepare a regular matrix multiplication.
883 assert p*q == 1
884 other = other*cvxopt.spdiag([1]*n)
885 p, q = other._shape
887 assert n == p
889 # Handle the common row by column multiplication more efficiently.
890 # This includes some scalar by scalar products (both sides nonconstant).
891 row_by_column = m == 1 and q == 1
893 # Regular matrix multiplication.
894 coefs = {}
895 for (lmtbs, lcoef), (rmtbs, rcoef) in product(
896 self._coefs.items(), other._coefs.items()):
897 if len(lmtbs) + len(rmtbs) > 2:
898 raise NotImplementedError("Product not supported: "
899 "One operand is biaffine and the other nonconstant.")
900 elif len(lmtbs) == 0:
901 # Constant by constant, linear or bilinear.
902 if row_by_column:
903 factor = lcoef.T
904 else:
905 # Use identity vec(AB) = (I ⊗ A)vec(B).
906 factor = left_kronecker_I(lcoef, q, reshape=(m, n))
908 self._update_coefs(coefs, rmtbs, factor*rcoef)
909 elif len(rmtbs) == 0:
910 # Constant, linear or bilinear by constant.
911 if row_by_column:
912 factor = rcoef.T
913 else:
914 # Use identity vec(AB) = (Bᵀ ⊗ I)vec(A).
915 factor = right_kronecker_I(
916 rcoef, m, reshape=(p, q), postT=True)
918 self._update_coefs(coefs, lmtbs, factor*lcoef)
919 else:
920 # Linear by linear.
921 assert len(lmtbs) == 1 and len(rmtbs) == 1
923 if row_by_column:
924 # Use identity (Ax)ᵀ(By) = vec(AᵀB)ᵀ·vec(xyᵀ).
925 coef = (lcoef.T*rcoef)[:].T
926 else:
927 # Recipe found in "Robust conic optimization in Python"
928 # (Stahlberg 2020).
929 a, b = lmtbs[0].dim, rmtbs[0].dim
931 A = cvxopt_K(m, n)*lcoef
932 B = rcoef
934 A = A.T
935 B = B.T
936 A.size = m*n*a, 1
937 B.size = p*q*b, 1
939 A = cvxopt.spdiag([cvxopt_K(a, n)]*m)*A
940 B = cvxopt.spdiag([cvxopt_K(b, p)]*q)*B
941 A.size = n, m*a
942 B.size = p, q*b
943 A = A.T
945 C = A*B
946 C.size = a, m*q*b
947 C = C*cvxopt.spdiag([cvxopt_K(b, m)]*q)
948 C.size = a*b, m*q
949 C = C.T
951 coef = C
953 # Forbid quadratic results.
954 if coef and lmtbs[0] is rmtbs[0]:
955 raise TypeError("The product of {} and {} is quadratic "
956 "in {} but a biaffine result is required here."
957 .format(self.string, other.string, lmtbs[0].name))
959 self._update_coefs(coefs, lmtbs + rmtbs, coef)
961 return self._common_basetype(other)(string, (m, q), coefs)
963 @convert_operands(lMatMul=True)
964 @refine_operands(stop_at_affine=True)
965 def __rmul__(self, other):
966 """Denote matrix multiplication from the left hand side."""
967 if not isinstance(other, BiaffineExpression):
968 return NotImplemented
970 return other.__mul__(self)
972 @convert_operands(sameShape=True)
973 @refine_operands(stop_at_affine=True)
974 def __or__(self, other):
975 r"""Denote the scalar product with self on the left hand side.
977 For (complex) vectors :math:`a` and :math:`b` this is the dot product
979 .. math::
980 (a \mid b)
981 &= \langle a, b \rangle \\
982 &= a \cdot b \\
983 &= b^H a.
985 For (complex) matrices :math:`A` and :math:`B` this is the Frobenius
986 inner product
988 .. math::
989 (A \mid B)
990 &= \langle A, B \rangle_F \\
991 &= A : B \\
992 &= \operatorname{tr}(B^H A) \\
993 &= \operatorname{vec}(B)^H \operatorname{vec}(\overline{A})
995 .. note::
996 Write ``(A|B)`` instead of ``A|B`` for the scalar product of ``A``
997 and ``B`` to obtain correct operator binding within a larger
998 expression context.
999 """
1000 if not isinstance(other, BiaffineExpression):
1001 return NotImplemented
1003 result = other.vec.H * self.vec
1004 result._symbStr = glyphs.clever_dotp(
1005 self.string, other.string, other.complex, self.scalar)
1006 return result
1008 @convert_operands(sameShape=True)
1009 @refine_operands(stop_at_affine=True)
1010 def __ror__(self, other):
1011 """Denote the scalar product with self on the right hand side."""
1012 if not isinstance(other, BiaffineExpression):
1013 return NotImplemented
1015 return other.__or__(self)
1017 @convert_operands(sameShape=True)
1018 @refine_operands(stop_at_affine=True)
1019 def __xor__(self, other):
1020 """Denote the elementwise (or Hadamard) product."""
1021 if not isinstance(other, BiaffineExpression):
1022 return NotImplemented
1024 string = glyphs.hadamard(self.string, other.string)
1026 if other.constant:
1027 factor = cvxopt.spdiag(other._constant_coef)
1028 coefs = {mtbs: factor*coef for mtbs, coef in self._coefs.items()}
1029 elif self.constant:
1030 factor = cvxopt.spdiag(self._constant_coef)
1031 coefs = {mtbs: factor*coef for mtbs, coef in other._coefs.items()}
1032 else:
1033 return (self.diag*other.vec).reshaped(self._shape).renamed(string)
1035 return self._common_basetype(other)(string, self._shape, coefs)
1037 @convert_operands(sameShape=True)
1038 @refine_operands(stop_at_affine=True)
1039 def __rxor__(self, other):
1040 """Denote the elementwise (or Hadamard) product."""
1041 if not isinstance(other, BiaffineExpression):
1042 return NotImplemented
1044 return other.__xor__(self)
1046 @convert_operands()
1047 @refine_operands(stop_at_affine=True)
1048 def __matmul__(self, other):
1049 """Denote the Kronecker product from the right hand side."""
1050 if not isinstance(other, BiaffineExpression):
1051 return NotImplemented
1053 cls = self._common_basetype(other)
1054 tc = cls._get_typecode()
1056 string = glyphs.kron(self.string, other.string)
1058 m, n = self._shape
1059 p, q = other._shape
1061 # Recipe in "Robust conic optimization in Python" (Stahlberg 2020).
1062 Kqm = cvxopt_K(q, m)
1063 Kqm_Ip = right_kronecker_I(Kqm, p)
1064 In_Kqm_Ip = left_kronecker_I(Kqm_Ip, n)
1066 def _kron(A, B):
1067 A_B = load_data(numpy.kron(cvx2np(A), cvx2np(B)), typecode=tc)[0]
1068 Kab = cvxopt_K(A.size[1], B.size[1])
1069 return In_Kqm_Ip*A_B*Kab
1071 coefs = {}
1072 for (lmtbs, lcoef), (rmtbs, rcoef) in product(
1073 self._coefs.items(), other._coefs.items()):
1074 if len(lmtbs) + len(rmtbs) <= 2:
1075 if len(lmtbs) == 1 and len(rmtbs) == 1 and lmtbs[0] is rmtbs[0]:
1076 raise TypeError("The product of {} and {} is quadratic "
1077 "in {} but a biaffine result is required here."
1078 .format(self.string, other.string, lmtbs[0].name))
1080 self._update_coefs(coefs, lmtbs + rmtbs, _kron(lcoef, rcoef))
1081 else:
1082 raise NotImplementedError("Product not supported: "
1083 "One operand is biaffine and the other nonconstant.")
1085 return cls(string, (m*p, n*q), coefs)
1087 @convert_operands()
1088 @refine_operands(stop_at_affine=True)
1089 def __rmatmul__(self, other):
1090 """Denote the Kronecker product from the left hand side."""
1091 if not isinstance(other, BiaffineExpression):
1092 return NotImplemented
1094 return other.__matmul__(self)
1096 @deprecated("2.2", "Use infix @ instead.")
1097 def kron(self, other):
1098 """Denote the Kronecker product from the right hand side."""
1099 return self.__matmul__(other)
1101 @deprecated("2.2", "Reverse operands and use infix @ instead.")
1102 def leftkron(self, other):
1103 """Denote the Kronecker product from the left hand side."""
1104 return self.__rmatmul__(other)
1106 @convert_operands(scalarRHS=True)
1107 @refine_operands(stop_at_affine=True)
1108 def __truediv__(self, other):
1109 """Denote elementwise division by a scalar constant."""
1110 if not isinstance(other, BiaffineExpression):
1111 return NotImplemented
1113 if not other.constant:
1114 raise TypeError(
1115 "You may only divide {} by a constant.".format(self.string))
1117 if other.is0:
1118 raise ZeroDivisionError(
1119 "Tried to divide {} by zero.".format(self.string))
1121 divisor = other._constant_coef[0]
1123 string = glyphs.div(self.string, other.string)
1124 coefs = {mtbs: coef / divisor for mtbs, coef in self._coefs.items()}
1126 return self._common_basetype(other)(string, self._shape, coefs)
1128 @convert_operands(scalarLHS=True)
1129 @refine_operands(stop_at_affine=True)
1130 def __rtruediv__(self, other):
1131 """Denote elementwise division by scalar constant self."""
1132 if not isinstance(other, BiaffineExpression):
1133 return NotImplemented
1135 return other.__truediv__(self)
1137 @convert_operands(horiCat=True)
1138 @refine_operands(stop_at_affine=True)
1139 def __and__(self, other):
1140 """Denote horizontal concatenation from the right hand side."""
1141 if not isinstance(other, BiaffineExpression):
1142 return NotImplemented
1144 string = glyphs.matrix_cat(self.string, other.string, horizontal=True)
1145 shape = (self._shape[0], self._shape[1] + other._shape[1])
1147 coefs = {}
1148 for mtbs in set(self._coefs.keys()).union(other._coefs.keys()):
1149 l = self._coefs[mtbs] if mtbs in self._coefs else cvxopt.spmatrix(
1150 [], [], [], (len(self), other._coefs[mtbs].size[1]))
1151 r = other._coefs[mtbs] if mtbs in other._coefs else cvxopt.spmatrix(
1152 [], [], [], (len(other), self._coefs[mtbs].size[1]))
1154 coefs[mtbs] = cvxopt_vcat([l, r])
1156 return self._common_basetype(other)(string, shape, coefs)
1158 @convert_operands(horiCat=True)
1159 @refine_operands(stop_at_affine=True)
1160 def __rand__(self, other):
1161 """Denote horizontal concatenation from the left hand side."""
1162 if not isinstance(other, BiaffineExpression):
1163 return NotImplemented
1165 return other.__and__(self)
1167 @convert_operands(vertCat=True)
1168 @refine_operands(stop_at_affine=True)
1169 def __floordiv__(self, other):
1170 """Denote vertical concatenation from below."""
1171 def interleave_columns(upper, lower, upperRows, lowerRows, cols):
1172 p, q = upperRows, lowerRows
1173 return [column for columnPairs in [
1174 (upper[j*p:j*p+p, :], lower[j*q:j*q+q, :]) for j in range(cols)]
1175 for column in columnPairs]
1177 if not isinstance(other, BiaffineExpression):
1178 return NotImplemented
1180 string = glyphs.matrix_cat(self.string, other.string, horizontal=False)
1182 p, q, c = self._shape[0], other._shape[0], self._shape[1]
1183 shape = (p + q, c)
1185 coefs = {}
1186 for mtbs in set(self._coefs.keys()).union(other._coefs.keys()):
1187 u = self._coefs[mtbs] if mtbs in self._coefs else cvxopt.spmatrix(
1188 [], [], [], (len(self), other._coefs[mtbs].size[1]))
1189 l = other._coefs[mtbs] if mtbs in other._coefs else cvxopt.spmatrix(
1190 [], [], [], (len(other), self._coefs[mtbs].size[1]))
1192 coefs[mtbs] = cvxopt_vcat(interleave_columns(u, l, p, q, c))
1194 return self._common_basetype(other)(string, shape, coefs)
1196 @convert_operands(vertCat=True)
1197 @refine_operands(stop_at_affine=True)
1198 def __rfloordiv__(self, other):
1199 """Denote vertical concatenation from above."""
1200 if not isinstance(other, BiaffineExpression):
1201 return NotImplemented
1203 return other.__floordiv__(self)
1205 @convert_operands(scalarRHS=True)
1206 @refine_operands() # Refine both sides to real if possible.
1207 def __pow__(self, other):
1208 """Denote exponentiation."""
1209 from .exp_powtrace import PowerTrace
1211 if not self.scalar:
1212 raise TypeError("May only exponentiate a scalar expression.")
1214 if not other.constant:
1215 raise TypeError("The exponent must be constant.")
1217 if other.complex:
1218 raise TypeError("The exponent must be real-valued.")
1220 exponent = other.value
1222 if exponent == 2:
1223 return (self | self) # Works for complex base.
1224 else:
1225 return PowerTrace(self, exponent) # Errors on complex base.
1227 # --------------------------------------------------------------------------
1228 # Properties and functions that describe the expression.
1229 # --------------------------------------------------------------------------
1231 @cached_property
1232 def hermitian(self): # noqa (D402 thinks this includes a signature)
1233 """Whether the expression is a hermitian (or symmetric) matrix.
1235 Uses :data:`~picos.settings.RELATIVE_HERMITIANNESS_TOLERANCE`.
1237 If PICOS rejects your near-hermitian (near-symmetric) expression as not
1238 hermitian (not symmetric), you can use :meth:`hermitianized` to correct
1239 larger numeric errors or the effects of noisy data.
1240 """
1241 return self.equals(
1242 self.H, relTol=settings.RELATIVE_HERMITIANNESS_TOLERANCE)
1244 @property
1245 def is0(self):
1246 """Whether this is a constant scalar, vector or matrix of all zeros."""
1247 return not self._coefs
1249 @cached_property
1250 def is1(self):
1251 """Whether this is a constant scalar or vector of all ones."""
1252 # Must be a scalar or vector.
1253 if self._shape[0] != 1 and self._shape[1] != 1:
1254 return False
1256 # Must be constant.
1257 if not self.constant:
1258 return False
1260 # Constant term must be all ones.
1261 return not self._constant_coef - 1
1263 @cached_property
1264 def isI(self):
1265 """Whether this is a constant identity matrix."""
1266 m, n = self._shape
1268 # Must be a square matrix.
1269 if m != n:
1270 return False
1272 # Must be constant.
1273 if not self.constant:
1274 return False
1276 # Constant term must be the identity.
1277 return not self._constant_coef - cvxopt.spdiag([1]*m)[:]
1279 @cached_property
1280 def complex(self):
1281 """Whether the expression can be complex-valued."""
1282 # NOTE: The typecode check works around a bug in CVXOPT prior to 1.2.4.
1283 return any(coef.typecode == "z" and coef.imag()
1284 for coef in self._coefs.values())
1286 @property
1287 def isreal(self):
1288 """Whether the expression is always real-valued."""
1289 return not self.complex
1291 @convert_operands()
1292 def equals(self, other, absTol=None, relTol=None):
1293 """Check mathematical equality with another (bi)affine expression.
1295 The precise type of both (bi)affine expressions may differ. In
1296 particular, a :class:`~.exp_affine.ComplexAffineExpression` with real
1297 coefficients and constant term can be equal to an
1298 :class:`~.exp_affine.AffineExpression`.
1300 If the operand is not already a PICOS expression, an attempt is made to
1301 load it as a constant affine expression. In this case, no reshaping or
1302 broadcasting is used to bring the constant term to the same shape as
1303 this expression. In particular,
1305 - ``0`` refers to a scalar zero (see also :meth:`is0`),
1306 - lists and tuples are treated as column vectors and
1307 - algebraic strings must specify a shape (see
1308 :func:`~.data.load_data`).
1310 :param other: Another PICOS expression or a constant numeric data value
1311 supported by :func:`~.data.load_data`.
1312 :param absTol: As long as all absolute differences between scalar
1313 entries of the coefficient matrices and the constant terms being
1314 compared does not exceed this bound, consider the expressions equal.
1315 :param relTol: As long as all absolute differences between scalar
1316 entries of the coefficient matrices and the constant terms being
1317 compared divided by the maximum absolute value found in either term
1318 does not exceed this bound, consider the expressions equal.
1320 :Example:
1322 >>> from picos import Constant
1323 >>> A = Constant("A", 0, (5,5))
1324 >>> repr(A)
1325 '<5×5 Real Constant: A>'
1326 >>> A.is0
1327 True
1328 >>> A.equals(0)
1329 False
1330 >>> A.equals("|0|(5,5)")
1331 True
1332 >>> repr(A*1j)
1333 '<5×5 Complex Constant: A·1j>'
1334 >>> A.equals(A*1j)
1335 True
1336 """
1337 if self is other:
1338 return True
1340 if not isinstance(other, BiaffineExpression):
1341 return False
1343 if self._shape != other._shape:
1344 return False
1346 assert not any((mtbs[1], mtbs[0]) in other._bilinear_coefs
1347 for mtbs in self._bilinear_coefs), \
1348 "{} and {} store bilinear terms in a different order." \
1349 .format(self.string, other.string)
1351 for mtbs in other._coefs:
1352 if mtbs not in self._coefs:
1353 return False
1355 for mtbs, coef in self._coefs.items():
1356 if mtbs not in other._coefs:
1357 return False
1359 if not cvxopt_equals(coef, other._coefs[mtbs], absTol, relTol):
1360 return False
1362 return True
1364 # --------------------------------------------------------------------------
1365 # Methods and properties that return modified copies.
1366 # --------------------------------------------------------------------------
1368 def renamed(self, string):
1369 """Return the expression with a modified string description."""
1370 return self._basetype(string, self._shape, self._coefs)
1372 def reshaped(self, shape, order="F"):
1373 r"""Return the expression reshaped in the given order.
1375 The default indexing order is column-major. Given an :math:`m \times n`
1376 matrix, reshaping in the default order is a constant time operation
1377 while reshaping in row-major order requires :math:`O(mn)` time. However,
1378 the latter allows you to be consistent with NumPy, which uses C-order (a
1379 generalization of row-major) by default.
1381 :param str order:
1382 The indexing order to use when reshaping. Must be either ``"F"`` for
1383 Fortran-order (column-major) or ``"C"`` for C-order (row-major).
1385 :Example:
1387 >>> from picos import Constant
1388 >>> C = Constant("C", range(6), (2, 3))
1389 >>> print(C)
1390 [ 0.00e+00 2.00e+00 4.00e+00]
1391 [ 1.00e+00 3.00e+00 5.00e+00]
1392 >>> print(C.reshaped((3, 2)))
1393 [ 0.00e+00 3.00e+00]
1394 [ 1.00e+00 4.00e+00]
1395 [ 2.00e+00 5.00e+00]
1396 >>> print(C.reshaped((3, 2), order="C"))
1397 [ 0.00e+00 2.00e+00]
1398 [ 4.00e+00 1.00e+00]
1399 [ 3.00e+00 5.00e+00]
1400 """
1401 if order not in "FC":
1402 raise ValueError("Order must be given as 'F' or 'C'.")
1404 shape = load_shape(shape, wildcards=True)
1406 if shape == self._shape:
1407 return self
1408 elif shape == (None, None):
1409 return self
1411 length = len(self)
1413 if shape[0] is None:
1414 shape = (length // shape[1], shape[1])
1415 elif shape[1] is None:
1416 shape = (shape[0], length // shape[0])
1418 if shape[0]*shape[1] != length:
1419 raise ValueError("Can only reshape to a matrix of same size.")
1421 if order == "F":
1422 coefs = self._coefs
1423 reshaped_glyph = glyphs.reshaped
1424 else:
1425 m, n = self._shape
1426 p, q = shape
1428 K_old = cvxopt_K(m, n, self._typecode)
1429 K_new = cvxopt_K(q, p, self._typecode)
1430 R = K_new * K_old
1432 coefs = {mtbs: R * coef for mtbs, coef in self._coefs.items()}
1433 reshaped_glyph = glyphs.reshaprm
1435 string = reshaped_glyph(self.string, glyphs.shape(shape))
1436 return self._basetype(string, shape, coefs)
1438 def broadcasted(self, shape):
1439 """Return the expression broadcasted to the given shape.
1441 :Example:
1443 >>> from picos import Constant
1444 >>> C = Constant("C", range(6), (2, 3))
1445 >>> print(C)
1446 [ 0.00e+00 2.00e+00 4.00e+00]
1447 [ 1.00e+00 3.00e+00 5.00e+00]
1448 >>> print(C.broadcasted((6, 6)))
1449 [ 0.00e+00 2.00e+00 4.00e+00 0.00e+00 2.00e+00 4.00e+00]
1450 [ 1.00e+00 3.00e+00 5.00e+00 1.00e+00 3.00e+00 5.00e+00]
1451 [ 0.00e+00 2.00e+00 4.00e+00 0.00e+00 2.00e+00 4.00e+00]
1452 [ 1.00e+00 3.00e+00 5.00e+00 1.00e+00 3.00e+00 5.00e+00]
1453 [ 0.00e+00 2.00e+00 4.00e+00 0.00e+00 2.00e+00 4.00e+00]
1454 [ 1.00e+00 3.00e+00 5.00e+00 1.00e+00 3.00e+00 5.00e+00]
1455 """
1456 shape = load_shape(shape, wildcards=True)
1457 shape = blend_shapes(shape, self._shape)
1459 if shape == self._shape:
1460 return self
1462 vdup = shape[0] // self._shape[0]
1463 hdup = shape[1] // self._shape[1]
1465 if (self._shape[0] * vdup, self._shape[1] * hdup) != shape:
1466 raise ValueError("Cannot broadcast from shape {} to {}."
1467 .format(glyphs.shape(self._shape), glyphs.shape(shape)))
1469 if self._shape == (1, 1):
1470 string = glyphs.matrix(self.string)
1471 return (self * cvxopt.matrix(1.0, shape)).renamed(string)
1473 string = glyphs.bcasted(self.string, glyphs.shape(shape))
1474 return (cvxopt.matrix(1.0, (vdup, hdup)) @ self).renamed(string)
1476 def reshaped_or_broadcasted(self, shape):
1477 """Return the expression :meth:`reshaped` or :meth:`broadcasted`.
1479 Unlike with :meth:`reshaped` and :meth:`broadcasted`, the target shape
1480 may not contain a wildcard character.
1482 If the wildcard-free target shape has the same number of elements as
1483 the current shape, then this is the same as :meth:`reshaped`, otherwise
1484 it is the same as :meth:`broadcasted`.
1485 """
1486 shape = load_shape(shape, wildcards=False)
1488 try:
1489 if shape[0]*shape[1] == len(self):
1490 return self.reshaped(shape)
1491 else:
1492 return self.broadcasted(shape)
1493 except ValueError:
1494 raise ValueError(
1495 "Cannot reshape or broadcast from shape {} to {}.".format(
1496 glyphs.shape(self._shape), glyphs.shape(shape))) from None
1498 @cached_property
1499 def hermitianized(self):
1500 r"""The expression projected onto the subspace of hermitian matrices.
1502 For a square (complex) affine expression :math:`A`, this is
1503 :math:`\frac{1}{2}(A + A^H)`.
1505 If the expression is not complex, then this is a projection onto the
1506 subspace of symmetric matrices.
1507 """
1508 if not self.square:
1509 raise TypeError("Cannot hermitianize non-square {}.".format(self))
1511 return (self + self.H)/2
1513 @cached_property
1514 def real(self):
1515 """Real part of the expression."""
1516 return self._basetype(glyphs.real(self.string), self._shape,
1517 {mtbs: coef.real() for mtbs, coef in self._coefs.items()})
1519 @cached_property
1520 def imag(self):
1521 """Imaginary part of the expression."""
1522 return self._basetype(glyphs.imag(self.string), self._shape,
1523 {mtbs: coef.imag() for mtbs, coef in self._coefs.items()})
1525 @cached_property
1526 def bilin(self):
1527 """Pure bilinear part of the expression."""
1528 return self._basetype(
1529 glyphs.blinpart(self._symbStr), self._shape, self._bilinear_coefs)
1531 @cached_property
1532 def lin(self):
1533 """Linear part of the expression."""
1534 return self._basetype(
1535 glyphs.linpart(self._symbStr), self._shape, self._linear_coefs)
1537 @cached_property
1538 def noncst(self):
1539 """Nonconstant part of the expression."""
1540 coefs = {mtbs: coefs for mtbs, coefs in self._coefs.items() if mtbs}
1542 return self._basetype(
1543 glyphs.ncstpart(self._symbStr), self._shape, coefs)
1545 @cached_property
1546 def cst(self):
1547 """Constant part of the expression."""
1548 coefs = {(): self._coefs[()]} if () in self._coefs else {}
1550 return self._basetype(glyphs.cstpart(self._symbStr), self._shape, coefs)
1552 @cached_selfinverse_property
1553 def T(self):
1554 """Matrix transpose."""
1555 if len(self) == 1:
1556 return self
1558 m, n = self._shape
1559 K = cvxopt_K(m, n, self._typecode)
1561 string = glyphs.transp(self.string)
1562 shape = (self._shape[1], self._shape[0])
1563 coefs = {mtbs: K * coef for mtbs, coef in self._coefs.items()}
1565 return self._basetype(string, shape, coefs)
1567 @cached_selfinverse_property
1568 def conj(self):
1569 """Complex conjugate."""
1570 string = glyphs.conj(self.string)
1571 coefs = {mtbs: coef.H.T for mtbs, coef in self._coefs.items()}
1573 return self._basetype(string, self._shape, coefs)
1575 @cached_selfinverse_property
1576 def H(self):
1577 """Conjugate (or Hermitian) transpose."""
1578 return self.T.conj.renamed(glyphs.htransp(self._symbStr))
1580 def _square_equal_subsystem_dims(self, diagLen):
1581 """Support :func:`partial_trace` and :func:`partial_transpose`."""
1582 m, n = self._shape
1583 k = math.log(m, diagLen)
1585 if m != n or int(k) != k:
1586 raise TypeError("The expression has shape {} so it cannot be "
1587 "decomposed into subsystems of shape {}.".format(
1588 glyphs.shape(self._shape), glyphs.shape((diagLen,)*2)))
1590 return ((diagLen,)*2,)*int(k)
1592 def partial_transpose(self, subsystems, dimensions=2):
1593 r"""Return the expression with selected subsystems transposed.
1595 If the expression can be written as
1596 :math:`A_0 \otimes \cdots \otimes A_{n-1}` for matrices
1597 :math:`A_0, \ldots, A_{n-1}` with shapes given in ``dimensions``, then
1598 this returns :math:`B_0 \otimes \cdots \otimes B_{n-1}` with
1599 :math:`B_i = A_i^T`, if ``i in subsystems`` (with :math:`i = -1` read as
1600 :math:`n-1`), and :math:`B_i = A_i`, otherwise.
1602 :param subsystems: A collection of or a single subystem number, indexed
1603 from zero, corresponding to subsystems that shall be transposed.
1604 The value :math:`-1` refers to the last subsystem.
1605 :type subsystems: int or tuple or list
1607 :param dimensions: Either an integer :math:`d` so that the subsystems
1608 are assumed to be all of shape :math:`d \times d`, or a sequence of
1609 subsystem shapes where an integer :math:`d` within the sequence is
1610 read as :math:`d \times d`. In any case, the elementwise product
1611 over all subsystem shapes must equal the expression's shape.
1612 :type dimensions: int or tuple or list
1614 :raises TypeError: If the subsystems do not match the expression.
1615 :raises IndexError: If the subsystem selection is invalid.
1617 :Example:
1619 >>> from picos import Constant
1620 >>> A = Constant("A", range(16), (4, 4))
1621 >>> print(A) #doctest: +NORMALIZE_WHITESPACE
1622 [ 0.00e+00 4.00e+00 8.00e+00 1.20e+01]
1623 [ 1.00e+00 5.00e+00 9.00e+00 1.30e+01]
1624 [ 2.00e+00 6.00e+00 1.00e+01 1.40e+01]
1625 [ 3.00e+00 7.00e+00 1.10e+01 1.50e+01]
1626 >>> A0 = A.partial_transpose(0); A0
1627 <4×4 Real Constant: A.{[2×2]ᵀ⊗[2×2]}>
1628 >>> print(A0) #doctest: +NORMALIZE_WHITESPACE
1629 [ 0.00e+00 4.00e+00 2.00e+00 6.00e+00]
1630 [ 1.00e+00 5.00e+00 3.00e+00 7.00e+00]
1631 [ 8.00e+00 1.20e+01 1.00e+01 1.40e+01]
1632 [ 9.00e+00 1.30e+01 1.10e+01 1.50e+01]
1633 >>> A1 = A.partial_transpose(1); A1
1634 <4×4 Real Constant: A.{[2×2]⊗[2×2]ᵀ}>
1635 >>> print(A1) #doctest: +NORMALIZE_WHITESPACE
1636 [ 0.00e+00 1.00e+00 8.00e+00 9.00e+00]
1637 [ 4.00e+00 5.00e+00 1.20e+01 1.30e+01]
1638 [ 2.00e+00 3.00e+00 1.00e+01 1.10e+01]
1639 [ 6.00e+00 7.00e+00 1.40e+01 1.50e+01]
1640 """
1641 m, n = self._shape
1643 if isinstance(dimensions, int):
1644 dimensions = self._square_equal_subsystem_dims(dimensions)
1645 else:
1646 dimensions = [
1647 (d, d) if isinstance(d, int) else d for d in dimensions]
1649 if reduce(
1650 lambda x, y: (x[0]*y[0], x[1]*y[1]), dimensions) != self._shape:
1651 raise TypeError("Subsystem dimensions do not match expression.")
1653 if isinstance(subsystems, int):
1654 subsystems = (subsystems,)
1655 elif not subsystems:
1656 return self
1658 numSys = len(dimensions)
1659 subsystems = set(numSys - 1 if sys == -1 else sys for sys in subsystems)
1661 for sys in subsystems:
1662 if not isinstance(sys, int):
1663 raise IndexError("Subsystem indices must be integer, not {}."
1664 .format(type(sys).__name__))
1665 elif sys < 0:
1666 raise IndexError("Subsystem indices must be nonnegative.")
1667 elif sys >= numSys:
1668 raise IndexError(
1669 "Subsystem index {} out of range for {} systems total."
1670 .format(sys, numSys))
1672 # If all subsystems are transposed, this is regular transposition.
1673 if len(subsystems) == numSys:
1674 return self.T
1676 # Prepare sparse K such that K·vec(A) = vec(partial_transpose(A)).
1677 d = m * n
1678 V = [1]*d
1679 I = range(d)
1680 J = cvxopt.matrix(I)
1681 T = cvxopt.matrix(0, J.size)
1682 obh, obw = 1, 1
1683 sysStrings = None
1685 # Apply transpositions starting with the rightmost Kronecker factor.
1686 for sys in range(numSys - 1, -1, -1):
1687 # Shape of current system.
1688 p, q = dimensions[sys]
1689 sysString = glyphs.matrix(glyphs.shape((p, q)))
1691 # Height/width of "inner" blocks being moved, initially scalars.
1692 ibh, ibw = obh, obw
1694 # Heigh/width of "outer" blocks whose relative position is
1695 # maintained but that are subject to transposition independently.
1696 # In the last iteration this is the shape of the resulting matrix.
1697 obh *= p
1698 obw *= q
1700 # Only transpose selected subsystems.
1701 if sys not in subsystems:
1702 sysStrings = glyphs.kron(sysString, sysStrings) \
1703 if sysStrings else sysString
1704 continue
1705 else:
1706 sysStrings = glyphs.kron(glyphs.transp(sysString), sysStrings) \
1707 if sysStrings else glyphs.transp(sysString)
1709 # Shape of outer blocks after transposition.
1710 obhT, obwT = obw // ibw * ibh, obh // ibh * ibw
1712 # Shape of full matrix after transposition.
1713 mT, nT = m // obh * obhT, n // obw * obwT
1715 for vi in I:
1716 # Full matrix column and row.
1717 c, r = divmod(vi, m)
1719 # Outer block vertical index and row within outer block,
1720 # outer block horizontal index and column within outer block.
1721 obi, obr = divmod(r, obh)
1722 obj, obc = divmod(c, obw)
1724 # Inner block vertical index and row within inner block,
1725 # inner block horizontal index and column within inner block.
1726 ibi, ibr = divmod(obr, ibh)
1727 ibj, ibc = divmod(obc, ibw)
1729 # (1) ibi*ibw + ibc is column within the transposed outer block;
1730 # adding obj*obwT yields the column in the transposed matrix.
1731 # (2) ibj*ibh + ibr is row within the transposed outer block;
1732 # adding obi*obhT yields the row in the transposed matrix.
1733 # (3) tvi is index within the vectorized transposed matrix.
1734 tvi = (obj*obwT + ibi*ibw + ibc)*mT \
1735 + (obi*obhT + ibj*ibh + ibr)
1737 # Prepare the transposition.
1738 T[tvi] = J[vi]
1740 # Apply the transposition.
1741 J, T = T, J
1742 m, n, obh, obw = mT, nT, obhT, obwT
1744 # Finalize the partial commutation matrix.
1745 K = cvxopt.spmatrix(V, I, J, (d, d), self._typecode)
1747 string = glyphs.ptransp_(self.string, sysStrings)
1748 shape = (m, n)
1749 coefs = {mtbs: K * coef for mtbs, coef in self._coefs.items()}
1751 return self._basetype(string, shape, coefs)
1753 @cached_property
1754 def T0(self):
1755 r"""Expression with the first :math:`2 \times 2` subsystem transposed.
1757 Only available for a :math:`2^k \times 2^k` matrix with all subsystems
1758 of shape :math:`2 \times 2`. Use :meth:`partial_transpose` otherwise.
1759 """
1760 return self.partial_transpose(subsystems=0)
1762 @cached_property
1763 def T1(self):
1764 r"""Expression with the second :math:`2 \times 2` subsystem transposed.
1766 Only available for a :math:`2^k \times 2^k` matrix with all subsystems
1767 of shape :math:`2 \times 2`. Use :meth:`partial_transpose` otherwise.
1768 """
1769 return self.partial_transpose(subsystems=1)
1771 @cached_property
1772 def T2(self):
1773 r"""Expression with the third :math:`2 \times 2` subsystem transposed.
1775 Only available for a :math:`2^k \times 2^k` matrix with all subsystems
1776 of shape :math:`2 \times 2`. Use :meth:`partial_transpose` otherwise.
1777 """
1778 return self.partial_transpose(subsystems=2)
1780 @cached_property
1781 def T3(self):
1782 r"""Expression with the fourth :math:`2 \times 2` subsystem transposed.
1784 Only available for a :math:`2^k \times 2^k` matrix with all subsystems
1785 of shape :math:`2 \times 2`. Use :meth:`partial_transpose` otherwise.
1786 """
1787 return self.partial_transpose(subsystems=3)
1789 @cached_property
1790 def Tl(self):
1791 r"""Expression with the last :math:`2 \times 2` subsystem transposed.
1793 Only available for a :math:`2^k \times 2^k` matrix with all subsystems
1794 of shape :math:`2 \times 2`. Use :meth:`partial_transpose` otherwise.
1795 """
1796 return self.partial_transpose(subsystems=-1)
1798 @staticmethod
1799 def _reindex_F(indices, source, destination):
1800 """Convert indices between different tensor shapes in Fortran-order."""
1801 new = []
1802 offset = 0
1803 factor = 1
1805 for index, base in zip(indices, source):
1806 offset += factor*index
1807 factor *= base
1809 for base in destination:
1810 offset, remainder = divmod(offset, base)
1811 new.append(remainder)
1813 return tuple(new)
1815 @staticmethod
1816 def _reindex_C(indices, source, destination):
1817 """Convert indices between different tensor shapes in C-order."""
1818 new = []
1819 offset = 0
1820 factor = 1
1822 for index, base in zip(reversed(indices), reversed(source)):
1823 offset += factor*index
1824 factor *= base
1826 for base in reversed(destination):
1827 offset, remainder = divmod(offset, base)
1828 new.insert(0, remainder)
1830 return tuple(new)
1832 def reshuffled(self, permutation="ikjl", dimensions=None, order="C"):
1833 """Return the reshuffled or realigned expression.
1835 This operation works directly on matrices. However, it is equivalent to
1836 the following sequence of operations:
1838 1. The matrix is reshaped to a tensor with the given ``dimensions`` and
1839 according to ``order``.
1840 2. The tensor's axes are permuted according to ``permutation``.
1841 3. The tensor is reshaped back to the shape of the original matrix
1842 according to ``order``.
1844 For comparison, the following function applies the same operation to a
1845 2D NumPy :class:`~numpy:numpy.ndarray`:
1847 .. code::
1849 def reshuffle_numpy(matrix, permutation, dimensions, order):
1850 P = "{} -> {}".format("".join(sorted(permutation)), permutation)
1851 reshuffled = numpy.reshape(matrix, dimensions, order)
1852 reshuffled = numpy.einsum(P, reshuffled)
1853 return numpy.reshape(reshuffled, matrix.shape, order)
1855 :param permutation:
1856 A sequence of comparable elements with length equal to the number of
1857 tensor dimensions. The sequence is compared to its ordered version
1858 and the resulting permutation pattern is used to permute the tensor
1859 indices. For instance, the string ``"ikjl"`` is compared to its
1860 sorted version ``"ijkl"`` and denotes that the second and third axis
1861 should be swapped.
1862 :type permutation:
1863 str or tuple or list
1865 :param dimensions:
1866 If this is an integer sequence, then it defines the dimensions of
1867 the tensor. If this is :obj:`None`, then the tensor is assumed to be
1868 hypercubic and the number of dimensions is inferred from the
1869 ``permutation`` argument.
1870 :type dimensions:
1871 None or tuple or list
1873 :param str order:
1874 The indexing order to use for the virtual reshaping. Must be either
1875 ``"F"`` for Fortran-order (generalization of column-major) or
1876 ``"C"`` for C-order (generalization of row-major). Note that PICOS
1877 usually reshapes in Fortran-order while NumPy defaults to C-order.
1879 :Example:
1881 >>> from picos import Constant
1882 >>> A = Constant("A", range(16), (4, 4))
1883 >>> print(A) #doctest: +NORMALIZE_WHITESPACE
1884 [ 0.00e+00 4.00e+00 8.00e+00 1.20e+01]
1885 [ 1.00e+00 5.00e+00 9.00e+00 1.30e+01]
1886 [ 2.00e+00 6.00e+00 1.00e+01 1.40e+01]
1887 [ 3.00e+00 7.00e+00 1.10e+01 1.50e+01]
1888 >>> R = A.reshuffled(); R
1889 <4×4 Real Constant: shuffled(A,ikjl,C)>
1890 >>> print(R) #doctest: +NORMALIZE_WHITESPACE
1891 [ 0.00e+00 4.00e+00 1.00e+00 5.00e+00]
1892 [ 8.00e+00 1.20e+01 9.00e+00 1.30e+01]
1893 [ 2.00e+00 6.00e+00 3.00e+00 7.00e+00]
1894 [ 1.00e+01 1.40e+01 1.10e+01 1.50e+01]
1895 >>> A.reshuffled("ji").equals(A.T) # Regular transposition.
1896 True
1897 >>> A.reshuffled("3214").equals(A.T0) # Partial transposition (1).
1898 True
1899 >>> A.reshuffled("1432").equals(A.T1) # Partial transposition (2).
1900 True
1901 """
1902 m, n = self._shape
1903 mn = m*n
1905 # Load the permutation.
1906 ordered = sorted(permutation)
1907 P = dict(enumerate(ordered.index(element) for element in permutation))
1909 if len(set(P.values())) < len(P):
1910 raise ValueError("The sequence defining the permutation appears to "
1911 "contain duplicate elements.")
1913 assert not set(P.keys()).symmetric_difference(set(P.values()))
1915 numDims = len(P)
1917 # Load the dimensions.
1918 guessDimensions = dimensions is None
1920 if guessDimensions:
1921 dimensions = (int(mn**(1.0 / numDims)),)*numDims
1922 else:
1923 if len(dimensions) != numDims:
1924 raise ValueError("The number of indices does not match the "
1925 "number of dimensions.")
1927 if reduce(int.__mul__, dimensions, 1) != mn:
1928 raise TypeError("The {} matrix {} cannot be reshaped to a {} "
1929 "tensor.".format(glyphs.shape(self.shape), self.string,
1930 "hypercubic order {}".format(numDims) if guessDimensions
1931 else glyphs.size("", "").join(str(d) for d in dimensions)))
1933 # Load the indexing order.
1934 if order not in "FC":
1935 raise ValueError("Order must be given as 'F' or 'C'.")
1937 reindex = self._reindex_F if order == "F" else self._reindex_C
1939 # Nothing to do for the neutral permutation.
1940 if all(key == val for key, val in P.items()):
1941 return self
1943 # Create a sparse mtrix R such that R·vec(A) = vec(reshuffle(A)).
1944 V, I, J = [1]*mn, [], range(mn)
1945 for i in range(mn):
1946 (k, j) = divmod(i, m) # (j, k) are column-major matrix indices.
1947 indices = reindex((j, k), (m, n), dimensions)
1948 newIndices = tuple(indices[P[d]] for d in range(numDims))
1949 newDimensions = tuple(dimensions[P[d]] for d in range(numDims))
1950 (j, k) = reindex(newIndices, newDimensions, (m, n))
1951 I.append(k*m + j)
1952 R = cvxopt.spmatrix(V, I, J, (mn, mn), self._typecode)
1954 # Create the string.
1955 strArgs = [self.string, str(permutation).replace(" ", ""), order]
1957 if not guessDimensions:
1958 strArgs.insert(2, str(dimensions).replace(" ", ""))
1960 string = glyphs.shuffled(",".join(strArgs))
1962 # Finalize the new expression.
1963 shape = (m, n)
1964 coefs = {mtbs: R * coef for mtbs, coef in self._coefs.items()}
1966 return self._basetype(string, shape, coefs)
1968 @cached_property
1969 def sum(self):
1970 """Sum over all scalar elements of the expression."""
1971 # NOTE: glyphs.clever_dotp detects this case and uses the sum glyph.
1972 # NOTE: 1 on the right hand side in case self is complex.
1973 return (self | 1)
1975 @cached_property
1976 def tr(self):
1977 """Trace of a square expression."""
1978 if not self.square:
1979 raise TypeError("Cannot compute the trace of non-square {}."
1980 .format(self.string))
1982 # NOTE: glyphs.clever_dotp detects this case and uses the trace glyph.
1983 # NOTE: "I" on the right hand side in case self is complex.
1984 return (self | "I")
1986 def partial_trace(self, subsystems, dimensions=2):
1987 r"""Return the partial trace over selected subsystems.
1989 If the expression can be written as
1990 :math:`A_0 \otimes \cdots \otimes A_{n-1}` for matrices
1991 :math:`A_0, \ldots, A_{n-1}` with shapes given in ``dimensions``, then
1992 this returns :math:`B_0 \otimes \cdots \otimes B_{n-1}` with
1993 :math:`B_i = \operatorname{tr}(A_i)`, if ``i in subsystems`` (with
1994 :math:`i = -1` read as :math:`n-1`), and :math:`B_i = A_i`, otherwise.
1996 :param subsystems: A collection of or a single subystem number, indexed
1997 from zero, corresponding to subsystems that shall be traced over.
1998 The value :math:`-1` refers to the last subsystem.
1999 :type subsystems: int or tuple or list
2001 :param dimensions: Either an integer :math:`d` so that the subsystems
2002 are assumed to be all of shape :math:`d \times d`, or a sequence of
2003 subsystem shapes where an integer :math:`d` within the sequence is
2004 read as :math:`d \times d`. In any case, the elementwise product
2005 over all subsystem shapes must equal the expression's shape.
2006 :type dimensions: int or tuple or list
2008 :raises TypeError: If the subsystems do not match the expression or if
2009 a non-square subsystem is to be traced over.
2010 :raises IndexError: If the subsystem selection is invalid in any other
2011 way.
2013 :Example:
2015 >>> from picos import Constant
2016 >>> A = Constant("A", range(16), (4, 4))
2017 >>> print(A) #doctest: +NORMALIZE_WHITESPACE
2018 [ 0.00e+00 4.00e+00 8.00e+00 1.20e+01]
2019 [ 1.00e+00 5.00e+00 9.00e+00 1.30e+01]
2020 [ 2.00e+00 6.00e+00 1.00e+01 1.40e+01]
2021 [ 3.00e+00 7.00e+00 1.10e+01 1.50e+01]
2022 >>> A0 = A.partial_trace(0); A0
2023 <2×2 Real Constant: A.{tr([2×2])⊗[2×2]}>
2024 >>> print(A0) #doctest: +NORMALIZE_WHITESPACE
2025 [ 1.00e+01 1.80e+01]
2026 [ 1.20e+01 2.00e+01]
2027 >>> A1 = A.partial_trace(1); A1
2028 <2×2 Real Constant: A.{[2×2]⊗tr([2×2])}>
2029 >>> print(A1) #doctest: +NORMALIZE_WHITESPACE
2030 [ 5.00e+00 2.10e+01]
2031 [ 9.00e+00 2.50e+01]
2032 """
2033 # Shape of the original matrix.
2034 m, n = self._shape
2036 if isinstance(dimensions, int):
2037 dimensions = self._square_equal_subsystem_dims(dimensions)
2038 else:
2039 dimensions = [
2040 (d, d) if isinstance(d, int) else d for d in dimensions]
2042 if reduce(
2043 lambda x, y: (x[0]*y[0], x[1]*y[1]), dimensions) != self._shape:
2044 raise TypeError("Subsystem dimensions do not match expression.")
2046 if isinstance(subsystems, int):
2047 subsystems = (subsystems,)
2048 elif not subsystems:
2049 return self
2051 numSys = len(dimensions)
2052 subsystems = set(numSys - 1 if sys == -1 else sys for sys in subsystems)
2054 for sys in subsystems:
2055 if not isinstance(sys, int):
2056 raise IndexError("Subsystem indices must be integer, not {}."
2057 .format(type(sys).__name__))
2058 elif sys < 0:
2059 raise IndexError("Subsystem indices must be nonnegative.")
2060 elif sys >= numSys:
2061 raise IndexError(
2062 "Subsystem index {} out of range for {} systems total."
2063 .format(sys, numSys))
2064 elif dimensions[sys][0] != dimensions[sys][1]:
2065 raise TypeError(
2066 "Subsystem index {} refers to a non-square subsystem that "
2067 "cannot be traced over.".format(sys))
2069 # If all subsystems are traced over, this is the regular trace.
2070 if len(subsystems) == numSys:
2071 return self.tr
2073 # Prepare sparse T such that T·vec(A) = vec(partial_trace(A)).
2074 T = []
2076 # Compute factors of T, one for each subsystem being traced over.
2077 ibh, ibw = m, n
2078 sysStrings = None
2079 for sys in range(numSys):
2080 # Shape of current system.
2081 p, q = dimensions[sys]
2082 sysString = glyphs.matrix(glyphs.shape((p, q)))
2084 # Heigh/width of "outer" blocks whose relative position is
2085 # maintained but that are reduced independently to the size of an
2086 # inner block if the current system is to be traced over.
2087 obh, obw = ibh, ibw
2089 # Height/width of "inner" blocks that are summed if they are
2090 # main-diagonal blocks of an outer block.
2091 ibh, ibw = obh // p, obw // q
2093 # Only trace over selected subsystems.
2094 if sys not in subsystems:
2095 sysStrings = glyphs.kron(sysStrings, sysString) \
2096 if sysStrings else sysString
2097 continue
2098 else:
2099 sysStrings = glyphs.kron(sysStrings, glyphs.trace(sysString)) \
2100 if sysStrings else glyphs.trace(sysString)
2102 # Shape of new matrix.
2103 assert p == q
2104 mN, nN = m // p, n // p
2106 # Prepare one factor of T.
2107 oldLen = m * n
2108 newLen = mN * nN
2109 V, I, J = [1]*(newLen*p), [], []
2110 shape = (newLen, oldLen)
2112 # Determine the summands that make up each entry of the new matrix.
2113 for viN in range(newLen):
2114 # A column/row pair that represents a scalar in the new matrix
2115 # and a number of p scalars within different on-diagonal inner
2116 # blocks but within the same outer block in the old matrix.
2117 cN, rN = divmod(viN, mN)
2119 # Index pair (obi, obj) for outer block in question, row/column
2120 # pair (ibr, ibc) identifies the scalar within each inner block.
2121 obi, ibr = divmod(rN, ibh)
2122 obj, ibc = divmod(cN, ibw)
2124 # Collect summands for the current entry of the new matrix; one
2125 # scalar per on-diagonal inner block.
2126 for k in range(p):
2127 rO = obi*obh + k*ibh + ibr
2128 cO = obj*obw + k*ibw + ibc
2129 I.append(viN)
2130 J.append(cO*m + rO)
2132 # Store the operator that performs the current subsystem trace.
2133 T.insert(0, cvxopt.spmatrix(V, I, J, shape, self._typecode))
2135 # Make next iteration work on the new matrix.
2136 m, n = mN, nN
2138 # Multiply out the linear partial trace operator T.
2139 # TODO: Fast matrix multiplication dynamic program useful here?
2140 T = reduce(lambda A, B: A*B, T)
2142 string = glyphs.ptrace_(self.string, sysStrings)
2143 shape = (m, n)
2144 coefs = {mtbs: T * coef for mtbs, coef in self._coefs.items()}
2146 return self._basetype(string, shape, coefs)
2148 @cached_property
2149 def tr0(self):
2150 r"""Expression with the first :math:`2 \times 2` subsystem traced out.
2152 Only available for a :math:`2^k \times 2^k` matrix with all subsystems
2153 of shape :math:`2 \times 2`. Use :meth:`partial_trace` otherwise.
2154 """
2155 return self.partial_trace(subsystems=0)
2157 @cached_property
2158 def tr1(self):
2159 r"""Expression with the second :math:`2 \times 2` subsystem traced out.
2161 Only available for a :math:`2^k \times 2^k` matrix with all subsystems
2162 of shape :math:`2 \times 2`. Use :meth:`partial_trace` otherwise.
2163 """
2164 return self.partial_trace(subsystems=1)
2166 @cached_property
2167 def tr2(self):
2168 r"""Expression with the third :math:`2 \times 2` subsystem traced out.
2170 Only available for a :math:`2^k \times 2^k` matrix with all subsystems
2171 of shape :math:`2 \times 2`. Use :meth:`partial_trace` otherwise.
2172 """
2173 return self.partial_trace(subsystems=2)
2175 @cached_property
2176 def tr3(self):
2177 r"""Expression with the fourth :math:`2 \times 2` subsystem traced out.
2179 Only available for a :math:`2^k \times 2^k` matrix with all subsystems
2180 of shape :math:`2 \times 2`. Use :meth:`partial_trace` otherwise.
2181 """
2182 return self.partial_trace(subsystems=3)
2184 @cached_property
2185 def trl(self):
2186 r"""Expression with the last :math:`2 \times 2` subsystem traced out.
2188 Only available for a :math:`2^k \times 2^k` matrix with all subsystems
2189 of shape :math:`2 \times 2`. Use :meth:`partial_trace` otherwise.
2190 """
2191 return self.partial_trace(subsystems=-1)
2193 @cached_property
2194 def vec(self):
2195 """Column-major vectorization of the expression as a column vector.
2197 .. note::
2198 Given an expression ``A``, ``A.vec`` and ``A[:]`` produce the same
2199 result (up to its string description) but ``A.vec`` is faster and
2200 its result is cached.
2202 :Example:
2204 >>> from picos import Constant
2205 >>> A = Constant("A", [[1, 2], [3, 4]])
2206 >>> A.vec.equals(A[:])
2207 True
2208 >>> A[:] is A[:]
2209 False
2210 >>> A.vec is A.vec
2211 True
2212 """
2213 if self._shape[1] == 1:
2214 return self
2215 else:
2216 return self._basetype(
2217 glyphs.vec(self.string), (len(self), 1), self._coefs)
2219 def dupvec(self, n):
2220 """Return a (repeated) column-major vectorization of the expression.
2222 :param int n: Number of times to duplicate the vectorization.
2224 :returns: A column vector.
2226 :Example:
2228 >>> from picos import Constant
2229 >>> A = Constant("A", [[1, 2], [3, 4]])
2230 >>> A.dupvec(1) is A.vec
2231 True
2232 >>> A.dupvec(3).equals(A.vec // A.vec // A.vec)
2233 True
2234 """
2235 if not isinstance(n, int):
2236 raise TypeError("Number of copies must be integer.")
2238 if n < 1:
2239 raise ValueError("Number of copies must be positive.")
2241 if n == 1:
2242 return self.vec
2243 else:
2244 string = glyphs.vec(glyphs.comma(self.string, n))
2245 shape = (len(self)*n, 1)
2246 coefs = {mtbs: cvxopt_vcat([coef]*n)
2247 for mtbs, coef in self._coefs.items()}
2249 return self._basetype(string, shape, coefs)
2251 @cached_property
2252 def trilvec(self):
2253 r"""Column-major vectorization of the lower triangular part.
2255 :returns:
2256 A column vector of all elements :math:`A_{ij}` that satisfy
2257 :math:`i \geq j`.
2259 .. note::
2261 If you want a row-major vectorization instead, write ``A.T.triuvec``
2262 instead of ``A.trilvec``.
2264 :Example:
2266 >>> from picos import Constant
2267 >>> A = Constant("A", [[1, 2], [3, 4], [5, 6]])
2268 >>> print(A)
2269 [ 1.00e+00 2.00e+00]
2270 [ 3.00e+00 4.00e+00]
2271 [ 5.00e+00 6.00e+00]
2272 >>> print(A.trilvec)
2273 [ 1.00e+00]
2274 [ 3.00e+00]
2275 [ 5.00e+00]
2276 [ 4.00e+00]
2277 [ 6.00e+00]
2278 """
2279 m, n = self._shape
2281 if n == 1: # Column vector or scalar.
2282 return self
2283 elif m == 1: # Row vector.
2284 return self[0]
2286 # Build a transformation D such that D·vec(A) = trilvec(A).
2287 rows = [j*m + i for j in range(n) for i in range(m) if i >= j]
2288 d = len(rows)
2289 D = cvxopt.spmatrix([1]*d, range(d), rows, (d, len(self)))
2291 string = glyphs.trilvec(self.string)
2292 shape = (d, 1)
2293 coefs = {mtbs: D*coef for mtbs, coef in self._coefs.items()}
2295 return self._basetype(string, shape, coefs)
2297 @cached_property
2298 def triuvec(self):
2299 r"""Column-major vectorization of the upper triangular part.
2301 :returns:
2302 A column vector of all elements :math:`A_{ij}` that satisfy
2303 :math:`i \leq j`.
2305 .. note::
2307 If you want a row-major vectorization instead, write ``A.T.trilvec``
2308 instead of ``A.triuvec``.
2310 :Example:
2312 >>> from picos import Constant
2313 >>> A = Constant("A", [[1, 2, 3], [4, 5, 6]])
2314 >>> print(A)
2315 [ 1.00e+00 2.00e+00 3.00e+00]
2316 [ 4.00e+00 5.00e+00 6.00e+00]
2317 >>> print(A.triuvec)
2318 [ 1.00e+00]
2319 [ 2.00e+00]
2320 [ 5.00e+00]
2321 [ 3.00e+00]
2322 [ 6.00e+00]
2323 """
2324 m, n = self._shape
2326 if m == 1: # Row vector or scalar.
2327 return self
2328 elif n == 1: # Column vector.
2329 return self[0]
2331 # Build a transformation D such that D·vec(A) = triuvec(A).
2332 rows = [j*m + i for j in range(n) for i in range(m) if i <= j]
2333 d = len(rows)
2334 D = cvxopt.spmatrix([1]*d, range(d), rows, (d, len(self)))
2336 string = glyphs.triuvec(self.string)
2337 shape = (d, 1)
2338 coefs = {mtbs: D*coef for mtbs, coef in self._coefs.items()}
2340 return self._basetype(string, shape, coefs)
2342 @cached_property
2343 def svec(self):
2344 """An isometric vectorization of a symmetric or hermitian expression.
2346 In the real symmetric case
2348 - the vectorization format is precisely the one define in [svec]_,
2349 - the vectorization is isometric and isomorphic, and
2350 - this is the same vectorization as used internally by the
2351 :class:`~.variables.SymmetricVariable` class.
2353 In the complex hermitian case
2355 - the same format is used, now resulting in a complex vector,
2356 - the vectorization is isometric but **not** isomorphic as there are
2357 guaranteed zeros in the imaginary part of the vector, and
2358 - this is **not** the same vectorization as the isomorphic,
2359 real-valued one used by :class:`~.variables.HermitianVariable`.
2361 The reverse operation is denoted by :attr:`desvec` in either case.
2363 :raises TypeError:
2364 If the expression is not square.
2366 :raises ValueError:
2367 If the expression is not hermitian.
2368 """
2369 if not self.square:
2370 raise TypeError("Can only compute svec for a square matrix, not for"
2371 " the {} expression {}.".format(self._shape, self.string))
2372 elif not self.hermitian:
2373 raise ValueError("Cannot compute svec for the non-hermitian "
2374 "expression {}.".format(self.string))
2376 vec = SymmetricVectorization(self._shape)
2377 V = vec._full2special
2379 string = glyphs.svec(self.string)
2381 if self.isreal:
2382 vec = SymmetricVectorization(self._shape)
2383 V = vec._full2special
2384 coefs = {mtbs: V*coef for mtbs, coef in self._coefs.items()}
2385 result = self._basetype(string, (vec.dim, 1), coefs)
2386 else:
2387 # NOTE: We need to reproduce svec using triuvec because, for numeric
2388 # reasons, SymmetricVectorization averages the elements in the
2389 # lower and upper triangular part. For symmetric matrices,
2390 # this is equivalent to the formal definition of svec found in
2391 # Datorro's book. For hermitian matrices however it is not.
2392 real_svec = self.real.svec
2393 imag_svec = 2**0.5 * 1j * self.imag.triuvec
2395 result = (real_svec + imag_svec).renamed(string)
2397 with unlocked_cached_properties(result):
2398 result.desvec = self
2400 return result
2402 @cached_property
2403 def desvec(self):
2404 r"""The reverse operation of :attr:`svec`.
2406 :raises TypeError:
2407 If the expression is not a vector or has a dimension outside the
2408 permissible set
2410 .. math::
2412 \left\{ \frac{n(n + 1)}{2}
2413 \mid n \in \mathbb{Z}_{\geq 1} \right\}
2414 = \left\{ n \in \mathbb{Z}_{\geq 1}
2415 \mid \frac{1}{2} \left( \sqrt{8n + 1} - 1 \right)
2416 \in \mathbb{Z}_{\geq 1} \right\}.
2418 :raises ValueError:
2419 In the case of a complex vector, If the vector is not in the image
2420 of :attr:`svec`.
2421 """
2422 if 1 not in self.shape:
2423 raise TypeError(
2424 "Can only compute desvec for a vector, not for the {} "
2425 "expression {}.".format(glyphs.shape(self._shape), self.string))
2427 n = 0.5*((8*len(self) + 1)**0.5 - 1)
2428 if int(n) != n:
2429 raise TypeError("Cannot compute desvec for the {}-dimensional "
2430 "vector {} as this size is not a possible outcome of svec."
2431 .format(len(self), self.string))
2432 n = int(n)
2434 vec = SymmetricVectorization((n, n))
2435 D = vec._special2full
2436 string = glyphs.desvec(self.string)
2438 if self.isreal:
2439 coefs = {mtbs: D*coef for mtbs, coef in self._coefs.items()}
2440 result = self._basetype(string, (n, n), coefs)
2442 assert result.hermitian
2443 else:
2444 # NOTE: While :attr:`svec` performs essentially the same operation
2445 # for both symmetric and hermitian matrices, we now need to
2446 # treat the imaginary separately as the imaginary part of a
2447 # hermitian matrix is skew-symmetric instead of symmetric.
2448 V, I, J = [], D.I, D.J
2449 for v, i, j in zip(D.V, I, J):
2450 V.append(1j*v if i % n <= i // n else -1j*v)
2451 C = cvxopt.spmatrix(V, I, J, D.size, tc="z")
2453 real_desvec = self.real.desvec
2454 imag_desvec = self._basetype(string, (n, n), {
2455 mtbs: C*coef for mtbs, coef in self.imag._coefs.items()})
2457 result = (real_desvec + imag_desvec).renamed(string)
2459 if not result.hermitian:
2460 raise ValueError("The complex vector {} is not in the image of "
2461 "svec. Note that svec is not bijective in the complex case."
2462 .format(self.string))
2464 with unlocked_cached_properties(result):
2465 result.svec = self
2467 return result
2469 def dupdiag(self, n):
2470 """Return a matrix with the (repeated) expression on the diagonal.
2472 Vectorization is performed in column-major order.
2474 :param int n: Number of times to duplicate the vectorization.
2475 """
2476 # Vectorize and duplicate the expression.
2477 vec = self.dupvec(n)
2478 d = len(vec)
2480 # Build a transformation D such that D·vec(A) = vec(diag(vec(A))).
2481 ones = [1]*d
2482 D = cvxopt.spdiag(ones)[:]
2483 D = cvxopt.spmatrix(ones, D.I, range(d), (D.size[0], d))
2485 string = glyphs.diag(vec.string)
2486 shape = (d, d)
2487 coefs = {mtbs: D*coef for mtbs, coef in vec._coefs.items()}
2489 return self._basetype(string, shape, coefs)
2491 @cached_property
2492 def diag(self):
2493 """Diagonal matrix with the expression on the main diagonal.
2495 Vectorization is performed in column-major order.
2496 """
2497 return self.dupdiag(1)
2499 @cached_property
2500 def maindiag(self):
2501 """The main diagonal of the expression as a column vector."""
2502 if 1 in self._shape:
2503 return self[0]
2505 # Build a transformation D such that D·vec(A) = diag(A).
2506 step = self._shape[0] + 1
2507 rows = [i*step for i in range(min(self._shape))]
2508 d = len(rows)
2509 D = cvxopt.spmatrix([1]*d, range(d), rows, (d, len(self)))
2511 string = glyphs.maindiag(self.string)
2512 shape = (d, 1)
2513 coefs = {mtbs: D*coef for mtbs, coef in self._coefs.items()}
2515 return self._basetype(string, shape, coefs)
2517 def factor_out(self, mutable):
2518 r"""Factor out a single mutable from a vector expression.
2520 If this expression is a column vector :math:`a` that depends on some
2521 mutable :math:`x` with a trivial internal vectorization format (i.e.
2522 :class:`~.vectorizations.FullVectorization`), then this method, called
2523 with :math:`x` as its argument, returns a pair of expressions
2524 :math:`(a_x, a_0)` such that :math:`a = a_x\operatorname{vec}(x) + a_0`.
2526 :returns:
2527 Two refined :class:`BiaffineExpression` instances that do not depend
2528 on ``mutable``.
2530 :raises TypeError:
2531 If the expression is not a column vector or if ``mutable`` is not a
2532 :class:`~.mutable.Mutable` or does not have a trivial vectorization
2533 format.
2535 :raises LookupError:
2536 If the expression does not depend on ``mutable``.
2538 :Example:
2540 >>> from picos import RealVariable
2541 >>> from picos.uncertain import UnitBallPerturbationSet
2542 >>> x = RealVariable("x", 3)
2543 >>> z = UnitBallPerturbationSet("z", 3).parameter
2544 >>> a = ((2*x + 3)^z) + 4*x + 5; a
2545 <3×1 Uncertain Affine Expression: (2·x + [3])⊙z + 4·x + [5]>
2546 >>> sorted(a.mutables, key=lambda mtb: mtb.name)
2547 [<3×1 Real Variable: x>, <3×1 Perturbation: z>]
2548 >>> az, a0 = a.factor_out(z)
2549 >>> az
2550 <3×3 Real Affine Expression: ((2·x + [3])⊙z + 4·x + [5])_z>
2551 >>> a0
2552 <3×1 Real Affine Expression: ((2·x + [3])⊙z + 4·x + [5])_0>
2553 >>> sorted(az.mutables.union(a0.mutables), key=lambda mtb: mtb.name)
2554 [<3×1 Real Variable: x>]
2555 >>> (az*z + a0).equals(a)
2556 True
2557 """
2558 if self._shape[1] != 1:
2559 raise TypeError(
2560 "Can only factor out mutables from column vectors but {} has "
2561 "a shape of {}.".format(self.string, glyphs.shape(self._shape)))
2563 mtb = mutable
2565 if not isinstance(mtb, Mutable):
2566 raise TypeError("Can only factor out mutables, not instances of {}."
2567 .format(type(mtb).__name__))
2569 if not isinstance(mtb._vec, FullVectorization):
2570 raise TypeError("Can only factor out mutables with a trivial "
2571 "vectorization format but {} uses {}."
2572 .format(mtb.name, type(mtb._vec).__name__))
2574 if mtb not in self.mutables:
2575 raise LookupError("Cannot factor out {} from {} as the latter does "
2576 "not depend on the former.".format(mtb.name, self.string))
2578 ax_string = glyphs.index(self.string, mtb.name)
2579 ax_shape = (self._shape[0], mtb.dim)
2581 # Recipe in "Robust conic optimization in Python" (Stahlberg 2020).
2582 ax_coefs = {}
2583 for mtbs, coef in self._coefs.items():
2584 if mtb not in set(mtbs):
2585 continue
2586 elif len(mtbs) == 1:
2587 assert mtbs[0] is mtb
2588 self._update_coefs(ax_coefs, (), coef[:])
2589 else:
2590 if mtbs[0] is mtb:
2591 other_mtb = mtbs[1]
2592 C = coef*cvxopt_K(other_mtb.dim, mtb.dim)
2593 else:
2594 assert mtbs[1] is mtb
2595 other_mtb = mtbs[0]
2596 C = coef
2598 d = other_mtb.dim
2599 C = cvxopt.sparse([C[:, i*d:(i+1)*d] for i in range(mtb.dim)])
2600 self._update_coefs(ax_coefs, (other_mtb,), C)
2602 ax = self._basetype(ax_string, ax_shape, ax_coefs)
2604 a0_string = glyphs.index(self.string, 0)
2605 a0_shape = self._shape
2606 a0_coefs = {M: C for M, C in self._coefs.items() if mtb not in set(M)}
2607 a0 = self._basetype(a0_string, a0_shape, a0_coefs)
2609 return ax.refined, a0.refined
2611 # --------------------------------------------------------------------------
2612 # Backwards compatibility methods.
2613 # --------------------------------------------------------------------------
2615 @classmethod
2616 @deprecated("2.0", useInstead="from_constant")
2617 def fromScalar(cls, scalar):
2618 """Create a class instance from a numeric scalar."""
2619 return cls.from_constant(scalar, (1, 1))
2621 @classmethod
2622 @deprecated("2.0", useInstead="from_constant")
2623 def fromMatrix(cls, matrix, size=None):
2624 """Create a class instance from a numeric matrix."""
2625 return cls.from_constant(matrix, size)
2627 @deprecated("2.0", useInstead="object.__xor__")
2628 def hadamard(self, fact):
2629 """Denote the elementwise (or Hadamard) product."""
2630 return self.__xor__(fact)
2632 @deprecated("2.0", useInstead="~.expression.Expression.constant")
2633 def isconstant(self):
2634 """Whether the expression involves no mutables."""
2635 return self.constant
2637 @deprecated("2.0", useInstead="equals")
2638 def same_as(self, other):
2639 """Check mathematical equality with another affine expression."""
2640 return self.equals(other)
2642 @deprecated("2.0", useInstead="T")
2643 def transpose(self):
2644 """Return the matrix transpose."""
2645 return self.T
2647 @cached_property
2648 @deprecated("2.0", useInstead="partial_transpose", decoratorLevel=1)
2649 def Tx(self):
2650 """Auto-detect few subsystems of same shape and transpose the last."""
2651 m, n = self._shape
2652 dims = None
2654 for k in range(2, int(math.log(min(m, n), 2)) + 1):
2655 p, q = int(round(m**(1.0/k))), int(round(n**(1.0/k)))
2656 if m == p**k and n == q**k:
2657 dims = ((p, q),)*k
2658 break
2660 if dims:
2661 return self.partial_transpose(subsystems=-1, dimensions=dims)
2662 else:
2663 raise RuntimeError("Failed to auto-detect subsystem dimensions for "
2664 "partial transposition: Only supported for {} matrices, {}."
2665 .format(glyphs.shape(
2666 (glyphs.power("m", "k"), glyphs.power("n", "k"))),
2667 glyphs.ge("k", 2)))
2669 @deprecated("2.0", useInstead="conj")
2670 def conjugate(self):
2671 """Return the complex conjugate."""
2672 return self.conj
2674 @deprecated("2.0", useInstead="H")
2675 def Htranspose(self):
2676 """Return the conjugate (or Hermitian) transpose."""
2677 return self.H
2679 @deprecated("2.0", reason="PICOS expressions are now immutable.")
2680 def copy(self):
2681 """Return a deep copy of the expression."""
2682 from copy import copy as cp
2684 return self._basetype(cp(self._symbStr), self._shape,
2685 {mtbs: cp(coef) for mtbs, coef in self._coefs.items()})
2687 @deprecated("2.0", reason="PICOS expressions are now immutable.")
2688 def soft_copy(self):
2689 """Return a shallow copy of the expression."""
2690 return self._basetype(self._symbStr, self._shape, self._coefs)
2693# --------------------------------------
2694__all__ = api_end(_API_START, globals())