Coverage for picos/expressions/exp_biaffine.py: 82.57%
1205 statements
« prev ^ index » next coverage.py v6.5.0, created at 2023-02-15 14:21 +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) == 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) == 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 = other.vec.H * self.vec
970 result._symbStr = glyphs.clever_dotp(
971 self.string, other.string, other.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 A_B = load_data(numpy.kron(cvx2np(A), cvx2np(B)), typecode=tc)[0]
1030 Kab = cvxopt_K(A.size[1], B.size[1])
1031 return In_Kqm_Ip*A_B*Kab
1033 coefs = {}
1034 for (lmtbs, lcoef), (rmtbs, rcoef) in product(
1035 self._coefs.items(), other._coefs.items()):
1036 if len(lmtbs) + len(rmtbs) <= 2:
1037 if len(lmtbs) == 1 and len(rmtbs) == 1 and lmtbs[0] is rmtbs[0]:
1038 raise TypeError("The product of {} and {} is quadratic "
1039 "in {} but a biaffine result is required here."
1040 .format(self.string, other.string, lmtbs[0].name))
1042 self._update_coefs(coefs, lmtbs + rmtbs, _kron(lcoef, rcoef))
1043 else:
1044 raise NotImplementedError("Product not supported: "
1045 "One operand is biaffine and the other nonconstant.")
1047 return cls(string, (m*p, n*q), coefs)
1049 @convert_operands()
1050 @refine_operands(stop_at_affine=True)
1051 def __rmatmul__(self, other):
1052 if not isinstance(other, BiaffineExpression):
1053 return Expression.__rmatmul__(self, other)
1055 return other.__matmul__(self)
1057 @deprecated("2.2", "Use infix @ instead.")
1058 def kron(self, other):
1059 """Denote the Kronecker product with another expression on the right."""
1060 return self.__matmul__(other)
1062 @deprecated("2.2", "Reverse operands and use infix @ instead.")
1063 def leftkron(self, other):
1064 """Denote the Kronecker product with another expression on the left."""
1065 return self.__rmatmul__(other)
1067 @convert_operands(scalarRHS=True)
1068 @refine_operands(stop_at_affine=True)
1069 def __truediv__(self, other):
1070 if not isinstance(other, BiaffineExpression):
1071 return Expression.__truediv__(self, other)
1073 if not other.constant:
1074 raise TypeError(
1075 "You may only divide {} by a constant.".format(self.string))
1077 if other.is0:
1078 raise ZeroDivisionError(
1079 "Tried to divide {} by zero.".format(self.string))
1081 divisor = other._constant_coef[0]
1083 string = glyphs.div(self.string, other.string)
1084 coefs = {mtbs: coef / divisor for mtbs, coef in self._coefs.items()}
1086 return self._common_basetype(other)(string, self._shape, coefs)
1088 @convert_operands(scalarLHS=True)
1089 @refine_operands(stop_at_affine=True)
1090 def __rtruediv__(self, other):
1091 if not isinstance(other, BiaffineExpression):
1092 return Expression.__rtruediv__(self, other)
1094 return other.__truediv__(self)
1096 @convert_operands(horiCat=True)
1097 @refine_operands(stop_at_affine=True)
1098 def __and__(self, other):
1099 if not isinstance(other, BiaffineExpression):
1100 return Expression.__and__(self, other)
1102 string = glyphs.matrix_cat(self.string, other.string, horizontal=True)
1103 shape = (self._shape[0], self._shape[1] + other._shape[1])
1105 coefs = {}
1106 for mtbs in set(self._coefs.keys()).union(other._coefs.keys()):
1107 l = self._coefs[mtbs] if mtbs in self._coefs else cvxopt.spmatrix(
1108 [], [], [], (len(self), other._coefs[mtbs].size[1]))
1109 r = other._coefs[mtbs] if mtbs in other._coefs else cvxopt.spmatrix(
1110 [], [], [], (len(other), self._coefs[mtbs].size[1]))
1112 coefs[mtbs] = cvxopt_vcat([l, r])
1114 return self._common_basetype(other)(string, shape, coefs)
1116 @convert_operands(horiCat=True)
1117 @refine_operands(stop_at_affine=True)
1118 def __rand__(self, other):
1119 if not isinstance(other, BiaffineExpression):
1120 return Expression.__rand__(self, other)
1122 return other.__and__(self)
1124 @convert_operands(vertCat=True)
1125 @refine_operands(stop_at_affine=True)
1126 def __floordiv__(self, other):
1127 def interleave_columns(upper, lower, upperRows, lowerRows, cols):
1128 p, q = upperRows, lowerRows
1129 return [column for columnPairs in [
1130 (upper[j*p:j*p+p, :], lower[j*q:j*q+q, :]) for j in range(cols)]
1131 for column in columnPairs]
1133 if not isinstance(other, BiaffineExpression):
1134 return Expression.__floordiv__(self, other)
1136 string = glyphs.matrix_cat(self.string, other.string, horizontal=False)
1138 p, q, c = self._shape[0], other._shape[0], self._shape[1]
1139 shape = (p + q, c)
1141 coefs = {}
1142 for mtbs in set(self._coefs.keys()).union(other._coefs.keys()):
1143 u = self._coefs[mtbs] if mtbs in self._coefs else cvxopt.spmatrix(
1144 [], [], [], (len(self), other._coefs[mtbs].size[1]))
1145 l = other._coefs[mtbs] if mtbs in other._coefs else cvxopt.spmatrix(
1146 [], [], [], (len(other), self._coefs[mtbs].size[1]))
1148 coefs[mtbs] = cvxopt_vcat(interleave_columns(u, l, p, q, c))
1150 return self._common_basetype(other)(string, shape, coefs)
1152 @convert_operands(vertCat=True)
1153 @refine_operands(stop_at_affine=True)
1154 def __rfloordiv__(self, other):
1155 if not isinstance(other, BiaffineExpression):
1156 return Expression.__rfloordiv__(self, other)
1158 return other.__floordiv__(self)
1160 @convert_operands(scalarRHS=True)
1161 @refine_operands() # Refine both sides to real if possible.
1162 def __pow__(self, other):
1163 from .exp_powtrace import PowerTrace
1165 if not self.scalar:
1166 raise TypeError("May only exponentiate a scalar expression.")
1168 if not other.constant:
1169 raise TypeError("The exponent must be constant.")
1171 if other.complex:
1172 raise TypeError("The exponent must be real-valued.")
1174 exponent = other.value
1176 if exponent == 2:
1177 return (self | self) # Works for complex base.
1178 else:
1179 return PowerTrace(self, exponent) # Errors on complex base.
1181 # --------------------------------------------------------------------------
1182 # Properties and functions that describe the expression.
1183 # --------------------------------------------------------------------------
1185 @cached_property
1186 def hermitian(self): # noqa (D402 thinks this includes a signature)
1187 """Whether the expression is a hermitian (or symmetric) matrix.
1189 Uses :data:`~picos.settings.RELATIVE_HERMITIANNESS_TOLERANCE`.
1191 If PICOS rejects your near-hermitian (near-symmetric) expression as not
1192 hermitian (not symmetric), you can use :meth:`hermitianized` to correct
1193 larger numeric errors or the effects of noisy data.
1194 """
1195 return self.equals(
1196 self.H, relTol=settings.RELATIVE_HERMITIANNESS_TOLERANCE)
1198 @property
1199 def is0(self):
1200 """Whether this is a constant scalar, vector or matrix of all zeros."""
1201 return not self._coefs
1203 @cached_property
1204 def is1(self):
1205 """Whether this is a constant scalar or vector of all ones."""
1206 # Must be a scalar or vector.
1207 if self._shape[0] != 1 and self._shape[1] != 1:
1208 return False
1210 # Must be constant.
1211 if not self.constant:
1212 return False
1214 # Constant term must be all ones.
1215 return not self._constant_coef - 1
1217 @cached_property
1218 def isI(self):
1219 """Whether this is a constant identity matrix."""
1220 m, n = self._shape
1222 # Must be a square matrix.
1223 if m != n:
1224 return False
1226 # Must be constant.
1227 if not self.constant:
1228 return False
1230 # Constant term must be the identity.
1231 return not self._constant_coef - cvxopt.spdiag([1]*m)[:]
1233 @cached_property
1234 def complex(self):
1235 """Whether the expression can be complex-valued."""
1236 # NOTE: The typecode check works around a bug in CVXOPT prior to 1.2.4.
1237 return any(coef.typecode == "z" and coef.imag()
1238 for coef in self._coefs.values())
1240 @property
1241 def isreal(self):
1242 """Whether the expression is always real-valued."""
1243 return not self.complex
1245 @convert_operands()
1246 def equals(self, other, absTol=None, relTol=None):
1247 """Check mathematical equality with another (bi)affine expression.
1249 The precise type of both (bi)affine expressions may differ. In
1250 particular, a :class:`~.exp_affine.ComplexAffineExpression` with real
1251 coefficients and constant term can be equal to an
1252 :class:`~.exp_affine.AffineExpression`.
1254 If the operand is not already a PICOS expression, an attempt is made to
1255 load it as a constant affine expression. In this case, no reshaping or
1256 broadcasting is used to bring the constant term to the same shape as
1257 this expression. In particular,
1259 - ``0`` refers to a scalar zero (see also :meth:`is0`),
1260 - lists and tuples are treated as column vectors and
1261 - algebraic strings must specify a shape (see
1262 :func:`~.data.load_data`).
1264 :param other: Another PICOS expression or a constant numeric data value
1265 supported by :func:`~.data.load_data`.
1266 :param absTol: As long as all absolute differences between scalar
1267 entries of the coefficient matrices and the constant terms being
1268 compared does not exceed this bound, consider the expressions equal.
1269 :param relTol: As long as all absolute differences between scalar
1270 entries of the coefficient matrices and the constant terms being
1271 compared divided by the maximum absolute value found in either term
1272 does not exceed this bound, consider the expressions equal.
1274 :Example:
1276 >>> from picos import Constant
1277 >>> A = Constant("A", 0, (5,5))
1278 >>> repr(A)
1279 '<5×5 Real Constant: A>'
1280 >>> A.is0
1281 True
1282 >>> A.equals(0)
1283 False
1284 >>> A.equals("|0|(5,5)")
1285 True
1286 >>> repr(A*1j)
1287 '<5×5 Complex Constant: A·1j>'
1288 >>> A.equals(A*1j)
1289 True
1290 """
1291 if self is other:
1292 return True
1294 if not isinstance(other, BiaffineExpression):
1295 return False
1297 if self._shape != other._shape:
1298 return False
1300 assert not any((mtbs[1], mtbs[0]) in other._bilinear_coefs
1301 for mtbs in self._bilinear_coefs), \
1302 "{} and {} store bilinear terms in a different order." \
1303 .format(self.string, other.string)
1305 for mtbs in other._coefs:
1306 if mtbs not in self._coefs:
1307 return False
1309 for mtbs, coef in self._coefs.items():
1310 if mtbs not in other._coefs:
1311 return False
1313 if not cvxopt_equals(coef, other._coefs[mtbs], absTol, relTol):
1314 return False
1316 return True
1318 # --------------------------------------------------------------------------
1319 # Methods and properties that return modified copies.
1320 # --------------------------------------------------------------------------
1322 def renamed(self, string):
1323 """Return the expression with a modified string description."""
1324 return self._basetype(string, self._shape, self._coefs)
1326 def reshaped(self, shape, order="F"):
1327 r"""Return the expression reshaped in the given order.
1329 The default indexing order is column-major. Given an :math:`m \times n`
1330 matrix, reshaping in the default order is a constant time operation
1331 while reshaping in row-major order requires :math:`O(mn)` time. However,
1332 the latter allows you to be consistent with NumPy, which uses C-order (a
1333 generalization of row-major) by default.
1335 :param str order:
1336 The indexing order to use when reshaping. Must be either ``"F"`` for
1337 Fortran-order (column-major) or ``"C"`` for C-order (row-major).
1339 :Example:
1341 >>> from picos import Constant
1342 >>> C = Constant("C", range(6), (2, 3))
1343 >>> print(C)
1344 [ 0.00e+00 2.00e+00 4.00e+00]
1345 [ 1.00e+00 3.00e+00 5.00e+00]
1346 >>> print(C.reshaped((3, 2)))
1347 [ 0.00e+00 3.00e+00]
1348 [ 1.00e+00 4.00e+00]
1349 [ 2.00e+00 5.00e+00]
1350 >>> print(C.reshaped((3, 2), order="C"))
1351 [ 0.00e+00 2.00e+00]
1352 [ 4.00e+00 1.00e+00]
1353 [ 3.00e+00 5.00e+00]
1354 """
1355 if order not in "FC":
1356 raise ValueError("Order must be given as 'F' or 'C'.")
1358 shape = load_shape(shape, wildcards=True)
1360 if shape == self._shape:
1361 return self
1362 elif shape == (None, None):
1363 return self
1365 length = len(self)
1367 if shape[0] is None:
1368 shape = (length // shape[1], shape[1])
1369 elif shape[1] is None:
1370 shape = (shape[0], length // shape[0])
1372 if shape[0]*shape[1] != length:
1373 raise ValueError("Can only reshape to a matrix of same size.")
1375 if order == "F":
1376 coefs = self._coefs
1377 reshaped_glyph = glyphs.reshaped
1378 else:
1379 m, n = self._shape
1380 p, q = shape
1382 K_old = cvxopt_K(m, n, self._typecode)
1383 K_new = cvxopt_K(q, p, self._typecode)
1384 R = K_new * K_old
1386 coefs = {mtbs: R * coef for mtbs, coef in self._coefs.items()}
1387 reshaped_glyph = glyphs.reshaprm
1389 string = reshaped_glyph(self.string, glyphs.shape(shape))
1390 return self._basetype(string, shape, coefs)
1392 def broadcasted(self, shape):
1393 """Return the expression broadcasted to the given shape.
1395 :Example:
1397 >>> from picos import Constant
1398 >>> C = Constant("C", range(6), (2, 3))
1399 >>> print(C)
1400 [ 0.00e+00 2.00e+00 4.00e+00]
1401 [ 1.00e+00 3.00e+00 5.00e+00]
1402 >>> print(C.broadcasted((6, 6)))
1403 [ 0.00e+00 2.00e+00 4.00e+00 0.00e+00 2.00e+00 4.00e+00]
1404 [ 1.00e+00 3.00e+00 5.00e+00 1.00e+00 3.00e+00 5.00e+00]
1405 [ 0.00e+00 2.00e+00 4.00e+00 0.00e+00 2.00e+00 4.00e+00]
1406 [ 1.00e+00 3.00e+00 5.00e+00 1.00e+00 3.00e+00 5.00e+00]
1407 [ 0.00e+00 2.00e+00 4.00e+00 0.00e+00 2.00e+00 4.00e+00]
1408 [ 1.00e+00 3.00e+00 5.00e+00 1.00e+00 3.00e+00 5.00e+00]
1409 """
1410 shape = load_shape(shape, wildcards=True)
1411 shape = blend_shapes(shape, self._shape)
1413 if shape == self._shape:
1414 return self
1416 vdup = shape[0] // self._shape[0]
1417 hdup = shape[1] // self._shape[1]
1419 if (self._shape[0] * vdup, self._shape[1] * hdup) != shape:
1420 raise ValueError("Cannot broadcast from shape {} to {}."
1421 .format(glyphs.shape(self._shape), glyphs.shape(shape)))
1423 if self._shape == (1, 1):
1424 string = glyphs.matrix(self.string)
1425 return (self * cvxopt.matrix(1.0, shape)).renamed(string)
1427 string = glyphs.bcasted(self.string, glyphs.shape(shape))
1428 return (cvxopt.matrix(1.0, (vdup, hdup)) @ self).renamed(string)
1430 def reshaped_or_broadcasted(self, shape):
1431 """Return the expression :meth:`reshaped` or :meth:`broadcasted`.
1433 Unlike with :meth:`reshaped` and :meth:`broadcasted`, the target shape
1434 may not contain a wildcard character.
1436 If the wildcard-free target shape has the same number of elements as
1437 the current shape, then this is the same as :meth:`reshaped`, otherwise
1438 it is the same as :meth:`broadcasted`.
1439 """
1440 shape = load_shape(shape, wildcards=False)
1442 try:
1443 if shape[0]*shape[1] == len(self):
1444 return self.reshaped(shape)
1445 else:
1446 return self.broadcasted(shape)
1447 except ValueError:
1448 raise ValueError(
1449 "Cannot reshape or broadcast from shape {} to {}.".format(
1450 glyphs.shape(self._shape), glyphs.shape(shape))) from None
1452 @cached_property
1453 def hermitianized(self):
1454 r"""The expression projected onto the subspace of hermitian matrices.
1456 For a square (complex) affine expression :math:`A`, this is
1457 :math:`\frac{1}{2}(A + A^H)`.
1459 If the expression is not complex, then this is a projection onto the
1460 subspace of symmetric matrices.
1461 """
1462 if not self.square:
1463 raise TypeError("Cannot hermitianize non-square {}.".format(self))
1465 return (self + self.H)/2
1467 @cached_property
1468 def real(self):
1469 """Real part of the expression."""
1470 return self._basetype(glyphs.real(self.string), self._shape,
1471 {mtbs: coef.real() for mtbs, coef in self._coefs.items()})
1473 @cached_property
1474 def imag(self):
1475 """Imaginary part of the expression."""
1476 return self._basetype(glyphs.imag(self.string), self._shape,
1477 {mtbs: coef.imag() for mtbs, coef in self._coefs.items()})
1479 @cached_property
1480 def bilin(self):
1481 """Pure bilinear part of the expression."""
1482 return self._basetype(
1483 glyphs.blinpart(self._symbStr), self._shape, self._bilinear_coefs)
1485 @cached_property
1486 def lin(self):
1487 """Linear part of the expression."""
1488 return self._basetype(
1489 glyphs.linpart(self._symbStr), self._shape, self._linear_coefs)
1491 @cached_property
1492 def noncst(self):
1493 """Nonconstant part of the expression."""
1494 coefs = {mtbs: coefs for mtbs, coefs in self._coefs.items() if mtbs}
1496 return self._basetype(
1497 glyphs.ncstpart(self._symbStr), self._shape, coefs)
1499 @cached_property
1500 def cst(self):
1501 """Constant part of the expression."""
1502 coefs = {(): self._coefs[()]} if () in self._coefs else {}
1504 return self._basetype(glyphs.cstpart(self._symbStr), self._shape, coefs)
1506 @cached_selfinverse_property
1507 def T(self):
1508 """Matrix transpose."""
1509 if len(self) == 1:
1510 return self
1512 m, n = self._shape
1513 K = cvxopt_K(m, n, self._typecode)
1515 string = glyphs.transp(self.string)
1516 shape = (self._shape[1], self._shape[0])
1517 coefs = {mtbs: K * coef for mtbs, coef in self._coefs.items()}
1519 return self._basetype(string, shape, coefs)
1521 @cached_selfinverse_property
1522 def conj(self):
1523 """Complex conjugate."""
1524 string = glyphs.conj(self.string)
1525 coefs = {mtbs: coef.H.T for mtbs, coef in self._coefs.items()}
1527 return self._basetype(string, self._shape, coefs)
1529 @cached_selfinverse_property
1530 def H(self):
1531 """Conjugate (or Hermitian) transpose."""
1532 return self.T.conj.renamed(glyphs.htransp(self._symbStr))
1534 def _square_equal_subsystem_dims(self, diagLen):
1535 """Support :func:`partial_trace` and :func:`partial_transpose`."""
1536 m, n = self._shape
1537 k = math.log(m, diagLen)
1539 if m != n or int(k) != k:
1540 raise TypeError("The expression has shape {} so it cannot be "
1541 "decomposed into subsystems of shape {}.".format(
1542 glyphs.shape(self._shape), glyphs.shape((diagLen,)*2)))
1544 return ((diagLen,)*2,)*int(k)
1546 def partial_transpose(self, subsystems, dimensions=2):
1547 r"""Return the expression with selected subsystems transposed.
1549 If the expression can be written as
1550 :math:`A_0 \otimes \cdots \otimes A_{n-1}` for matrices
1551 :math:`A_0, \ldots, A_{n-1}` with shapes given in ``dimensions``, then
1552 this returns :math:`B_0 \otimes \cdots \otimes B_{n-1}` with
1553 :math:`B_i = A_i^T`, if ``i in subsystems`` (with :math:`i = -1` read as
1554 :math:`n-1`), and :math:`B_i = A_i`, otherwise.
1556 :param subsystems: A collection of or a single subystem number, indexed
1557 from zero, corresponding to subsystems that shall be transposed.
1558 The value :math:`-1` refers to the last subsystem.
1559 :type subsystems: int or tuple or list
1561 :param dimensions: Either an integer :math:`d` so that the subsystems
1562 are assumed to be all of shape :math:`d \times d`, or a sequence of
1563 subsystem shapes where an integer :math:`d` within the sequence is
1564 read as :math:`d \times d`. In any case, the elementwise product
1565 over all subsystem shapes must equal the expression's shape.
1566 :type dimensions: int or tuple or list
1568 :raises TypeError: If the subsystems do not match the expression.
1569 :raises IndexError: If the subsystem selection is invalid.
1571 :Example:
1573 >>> from picos import Constant
1574 >>> A = Constant("A", range(16), (4, 4))
1575 >>> print(A) #doctest: +NORMALIZE_WHITESPACE
1576 [ 0.00e+00 4.00e+00 8.00e+00 1.20e+01]
1577 [ 1.00e+00 5.00e+00 9.00e+00 1.30e+01]
1578 [ 2.00e+00 6.00e+00 1.00e+01 1.40e+01]
1579 [ 3.00e+00 7.00e+00 1.10e+01 1.50e+01]
1580 >>> A0 = A.partial_transpose(0); A0
1581 <4×4 Real Constant: A.{[2×2]ᵀ⊗[2×2]}>
1582 >>> print(A0) #doctest: +NORMALIZE_WHITESPACE
1583 [ 0.00e+00 4.00e+00 2.00e+00 6.00e+00]
1584 [ 1.00e+00 5.00e+00 3.00e+00 7.00e+00]
1585 [ 8.00e+00 1.20e+01 1.00e+01 1.40e+01]
1586 [ 9.00e+00 1.30e+01 1.10e+01 1.50e+01]
1587 >>> A1 = A.partial_transpose(1); A1
1588 <4×4 Real Constant: A.{[2×2]⊗[2×2]ᵀ}>
1589 >>> print(A1) #doctest: +NORMALIZE_WHITESPACE
1590 [ 0.00e+00 1.00e+00 8.00e+00 9.00e+00]
1591 [ 4.00e+00 5.00e+00 1.20e+01 1.30e+01]
1592 [ 2.00e+00 3.00e+00 1.00e+01 1.10e+01]
1593 [ 6.00e+00 7.00e+00 1.40e+01 1.50e+01]
1594 """
1595 m, n = self._shape
1597 if isinstance(dimensions, int):
1598 dimensions = self._square_equal_subsystem_dims(dimensions)
1599 else:
1600 dimensions = [
1601 (d, d) if isinstance(d, int) else d for d in dimensions]
1603 if reduce(
1604 lambda x, y: (x[0]*y[0], x[1]*y[1]), dimensions) != self._shape:
1605 raise TypeError("Subsystem dimensions do not match expression.")
1607 if isinstance(subsystems, int):
1608 subsystems = (subsystems,)
1609 elif not subsystems:
1610 return self
1612 numSys = len(dimensions)
1613 subsystems = set(numSys - 1 if sys == -1 else sys for sys in subsystems)
1615 for sys in subsystems:
1616 if not isinstance(sys, int):
1617 raise IndexError("Subsystem indices must be integer, not {}."
1618 .format(type(sys).__name__))
1619 elif sys < 0:
1620 raise IndexError("Subsystem indices must be nonnegative.")
1621 elif sys >= numSys:
1622 raise IndexError(
1623 "Subsystem index {} out of range for {} systems total."
1624 .format(sys, numSys))
1626 # If all subsystems are transposed, this is regular transposition.
1627 if len(subsystems) == numSys:
1628 return self.T
1630 # Prepare sparse K such that K·vec(A) = vec(partial_transpose(A)).
1631 d = m * n
1632 V = [1]*d
1633 I = range(d)
1634 J = cvxopt.matrix(I)
1635 T = cvxopt.matrix(0, J.size)
1636 obh, obw = 1, 1
1637 sysStrings = None
1639 # Apply transpositions starting with the rightmost Kronecker factor.
1640 for sys in range(numSys - 1, -1, -1):
1641 # Shape of current system.
1642 p, q = dimensions[sys]
1643 sysString = glyphs.matrix(glyphs.shape((p, q)))
1645 # Height/width of "inner" blocks being moved, initially scalars.
1646 ibh, ibw = obh, obw
1648 # Heigh/width of "outer" blocks whose relative position is
1649 # maintained but that are subject to transposition independently.
1650 # In the last iteration this is the shape of the resulting matrix.
1651 obh *= p
1652 obw *= q
1654 # Only transpose selected subsystems.
1655 if sys not in subsystems:
1656 sysStrings = glyphs.kron(sysString, sysStrings) \
1657 if sysStrings else sysString
1658 continue
1659 else:
1660 sysStrings = glyphs.kron(glyphs.transp(sysString), sysStrings) \
1661 if sysStrings else glyphs.transp(sysString)
1663 # Shape of outer blocks after transposition.
1664 obhT, obwT = obw // ibw * ibh, obh // ibh * ibw
1666 # Shape of full matrix after transposition.
1667 mT, nT = m // obh * obhT, n // obw * obwT
1669 for vi in I:
1670 # Full matrix column and row.
1671 c, r = divmod(vi, m)
1673 # Outer block vertical index and row within outer block,
1674 # outer block horizontal index and column within outer block.
1675 obi, obr = divmod(r, obh)
1676 obj, obc = divmod(c, obw)
1678 # Inner block vertical index and row within inner block,
1679 # inner block horizontal index and column within inner block.
1680 ibi, ibr = divmod(obr, ibh)
1681 ibj, ibc = divmod(obc, ibw)
1683 # (1) ibi*ibw + ibc is column within the transposed outer block;
1684 # adding obj*obwT yields the column in the transposed matrix.
1685 # (2) ibj*ibh + ibr is row within the transposed outer block;
1686 # adding obi*obhT yields the row in the transposed matrix.
1687 # (3) tvi is index within the vectorized transposed matrix.
1688 tvi = (obj*obwT + ibi*ibw + ibc)*mT \
1689 + (obi*obhT + ibj*ibh + ibr)
1691 # Prepare the transposition.
1692 T[tvi] = J[vi]
1694 # Apply the transposition.
1695 J, T = T, J
1696 m, n, obh, obw = mT, nT, obhT, obwT
1698 # Finalize the partial commutation matrix.
1699 K = cvxopt.spmatrix(V, I, J, (d, d), self._typecode)
1701 string = glyphs.ptransp_(self.string, sysStrings)
1702 shape = (m, n)
1703 coefs = {mtbs: K * coef for mtbs, coef in self._coefs.items()}
1705 return self._basetype(string, shape, coefs)
1707 @cached_property
1708 def T0(self):
1709 r"""Expression with the first :math:`2 \times 2` subsystem transposed.
1711 Only available for a :math:`2^k \times 2^k` matrix with all subsystems
1712 of shape :math:`2 \times 2`. Use :meth:`partial_transpose` otherwise.
1713 """
1714 return self.partial_transpose(subsystems=0)
1716 @cached_property
1717 def T1(self):
1718 r"""Expression with the second :math:`2 \times 2` subsystem transposed.
1720 Only available for a :math:`2^k \times 2^k` matrix with all subsystems
1721 of shape :math:`2 \times 2`. Use :meth:`partial_transpose` otherwise.
1722 """
1723 return self.partial_transpose(subsystems=1)
1725 @cached_property
1726 def T2(self):
1727 r"""Expression with the third :math:`2 \times 2` subsystem transposed.
1729 Only available for a :math:`2^k \times 2^k` matrix with all subsystems
1730 of shape :math:`2 \times 2`. Use :meth:`partial_transpose` otherwise.
1731 """
1732 return self.partial_transpose(subsystems=2)
1734 @cached_property
1735 def T3(self):
1736 r"""Expression with the fourth :math:`2 \times 2` subsystem transposed.
1738 Only available for a :math:`2^k \times 2^k` matrix with all subsystems
1739 of shape :math:`2 \times 2`. Use :meth:`partial_transpose` otherwise.
1740 """
1741 return self.partial_transpose(subsystems=3)
1743 @cached_property
1744 def Tl(self):
1745 r"""Expression with the last :math:`2 \times 2` subsystem transposed.
1747 Only available for a :math:`2^k \times 2^k` matrix with all subsystems
1748 of shape :math:`2 \times 2`. Use :meth:`partial_transpose` otherwise.
1749 """
1750 return self.partial_transpose(subsystems=-1)
1752 @staticmethod
1753 def _reindex_F(indices, source, destination):
1754 """Convert indices between different tensor shapes in Fortran-order."""
1755 new = []
1756 offset = 0
1757 factor = 1
1759 for index, base in zip(indices, source):
1760 offset += factor*index
1761 factor *= base
1763 for base in destination:
1764 offset, remainder = divmod(offset, base)
1765 new.append(remainder)
1767 return tuple(new)
1769 @staticmethod
1770 def _reindex_C(indices, source, destination):
1771 """Convert indices between different tensor shapes in C-order."""
1772 new = []
1773 offset = 0
1774 factor = 1
1776 for index, base in zip(reversed(indices), reversed(source)):
1777 offset += factor*index
1778 factor *= base
1780 for base in reversed(destination):
1781 offset, remainder = divmod(offset, base)
1782 new.insert(0, remainder)
1784 return tuple(new)
1786 def reshuffled(self, permutation="ikjl", dimensions=None, order="C"):
1787 """Return the reshuffled or realigned expression.
1789 This operation works directly on matrices. However, it is equivalent to
1790 the following sequence of operations:
1792 1. The matrix is reshaped to a tensor with the given ``dimensions`` and
1793 according to ``order``.
1794 2. The tensor's axes are permuted according to ``permutation``.
1795 3. The tensor is reshaped back to the shape of the original matrix
1796 according to ``order``.
1798 For comparison, the following function applies the same operation to a
1799 2D NumPy :class:`~numpy:numpy.ndarray`:
1801 .. code::
1803 def reshuffle_numpy(matrix, permutation, dimensions, order):
1804 P = "{} -> {}".format("".join(sorted(permutation)), permutation)
1805 reshuffled = numpy.reshape(matrix, dimensions, order)
1806 reshuffled = numpy.einsum(P, reshuffled)
1807 return numpy.reshape(reshuffled, matrix.shape, order)
1809 :param permutation:
1810 A sequence of comparable elements with length equal to the number of
1811 tensor dimensions. The sequence is compared to its ordered version
1812 and the resulting permutation pattern is used to permute the tensor
1813 indices. For instance, the string ``"ikjl"`` is compared to its
1814 sorted version ``"ijkl"`` and denotes that the second and third axis
1815 should be swapped.
1816 :type permutation:
1817 str or tuple or list
1819 :param dimensions:
1820 If this is an integer sequence, then it defines the dimensions of
1821 the tensor. If this is :obj:`None`, then the tensor is assumed to be
1822 hypercubic and the number of dimensions is inferred from the
1823 ``permutation`` argument.
1824 :type dimensions:
1825 None or tuple or list
1827 :param str order:
1828 The indexing order to use for the virtual reshaping. Must be either
1829 ``"F"`` for Fortran-order (generalization of column-major) or
1830 ``"C"`` for C-order (generalization of row-major). Note that PICOS
1831 usually reshapes in Fortran-order while NumPy defaults to C-order.
1833 :Example:
1835 >>> from picos import Constant
1836 >>> A = Constant("A", range(16), (4, 4))
1837 >>> print(A) #doctest: +NORMALIZE_WHITESPACE
1838 [ 0.00e+00 4.00e+00 8.00e+00 1.20e+01]
1839 [ 1.00e+00 5.00e+00 9.00e+00 1.30e+01]
1840 [ 2.00e+00 6.00e+00 1.00e+01 1.40e+01]
1841 [ 3.00e+00 7.00e+00 1.10e+01 1.50e+01]
1842 >>> R = A.reshuffled(); R
1843 <4×4 Real Constant: shuffled(A,ikjl,C)>
1844 >>> print(R) #doctest: +NORMALIZE_WHITESPACE
1845 [ 0.00e+00 4.00e+00 1.00e+00 5.00e+00]
1846 [ 8.00e+00 1.20e+01 9.00e+00 1.30e+01]
1847 [ 2.00e+00 6.00e+00 3.00e+00 7.00e+00]
1848 [ 1.00e+01 1.40e+01 1.10e+01 1.50e+01]
1849 >>> A.reshuffled("ji").equals(A.T) # Regular transposition.
1850 True
1851 >>> A.reshuffled("3214").equals(A.T0) # Partial transposition (1).
1852 True
1853 >>> A.reshuffled("1432").equals(A.T1) # Partial transposition (2).
1854 True
1855 """
1856 m, n = self._shape
1857 mn = m*n
1859 # Load the permutation.
1860 ordered = sorted(permutation)
1861 P = dict(enumerate(ordered.index(element) for element in permutation))
1863 if len(set(P.values())) < len(P):
1864 raise ValueError("The sequence defining the permutation appears to "
1865 "contain duplicate elements.")
1867 assert not set(P.keys()).symmetric_difference(set(P.values()))
1869 numDims = len(P)
1871 # Load the dimensions.
1872 guessDimensions = dimensions is None
1874 if guessDimensions:
1875 dimensions = (int(mn**(1.0 / numDims)),)*numDims
1876 else:
1877 if len(dimensions) != numDims:
1878 raise ValueError("The number of indices does not match the "
1879 "number of dimensions.")
1881 if reduce(int.__mul__, dimensions, 1) != mn:
1882 raise TypeError("The {} matrix {} cannot be reshaped to a {} "
1883 "tensor.".format(glyphs.shape(self.shape), self.string,
1884 "hypercubic order {}".format(numDims) if guessDimensions
1885 else glyphs.size("", "").join(str(d) for d in dimensions)))
1887 # Load the indexing order.
1888 if order not in "FC":
1889 raise ValueError("Order must be given as 'F' or 'C'.")
1891 reindex = self._reindex_F if order == "F" else self._reindex_C
1893 # Nothing to do for the neutral permutation.
1894 if all(key == val for key, val in P.items()):
1895 return self
1897 # Create a sparse mtrix R such that R·vec(A) = vec(reshuffle(A)).
1898 V, I, J = [1]*mn, [], range(mn)
1899 for i in range(mn):
1900 (k, j) = divmod(i, m) # (j, k) are column-major matrix indices.
1901 indices = reindex((j, k), (m, n), dimensions)
1902 newIndices = tuple(indices[P[d]] for d in range(numDims))
1903 newDimensions = tuple(dimensions[P[d]] for d in range(numDims))
1904 (j, k) = reindex(newIndices, newDimensions, (m, n))
1905 I.append(k*m + j)
1906 R = cvxopt.spmatrix(V, I, J, (mn, mn), self._typecode)
1908 # Create the string.
1909 strArgs = [self.string, str(permutation).replace(" ", ""), order]
1911 if not guessDimensions:
1912 strArgs.insert(2, str(dimensions).replace(" ", ""))
1914 string = glyphs.shuffled(",".join(strArgs))
1916 # Finalize the new expression.
1917 shape = (m, n)
1918 coefs = {mtbs: R * coef for mtbs, coef in self._coefs.items()}
1920 return self._basetype(string, shape, coefs)
1922 @cached_property
1923 def sum(self):
1924 """Sum over all scalar elements of the expression."""
1925 # NOTE: glyphs.clever_dotp detects this case and uses the sum glyph.
1926 # NOTE: 1 on the right hand side in case self is complex.
1927 return (self | 1)
1929 @cached_property
1930 def rowsum(self):
1931 """Sum over the rows of the expression as a row vector."""
1932 from .algebra import J
1934 return J(self.shape[0]).T * self
1936 @cached_property
1937 def colsum(self):
1938 """Sum over the columns of the expression as a column vector."""
1939 from .algebra import J
1941 return self * J(self.shape[1])
1943 @cached_property
1944 def tr(self):
1945 """Trace of a square expression."""
1946 if not self.square:
1947 raise TypeError("Cannot compute the trace of non-square {}."
1948 .format(self.string))
1950 # NOTE: glyphs.clever_dotp detects this case and uses the trace glyph.
1951 # NOTE: "I" on the right hand side in case self is complex.
1952 return (self | "I")
1954 def partial_trace(self, subsystems, dimensions=2):
1955 r"""Return the partial trace over selected subsystems.
1957 If the expression can be written as
1958 :math:`A_0 \otimes \cdots \otimes A_{n-1}` for matrices
1959 :math:`A_0, \ldots, A_{n-1}` with shapes given in ``dimensions``, then
1960 this returns :math:`B_0 \otimes \cdots \otimes B_{n-1}` with
1961 :math:`B_i = \operatorname{tr}(A_i)`, if ``i in subsystems`` (with
1962 :math:`i = -1` read as :math:`n-1`), and :math:`B_i = A_i`, otherwise.
1964 :param subsystems: A collection of or a single subystem number, indexed
1965 from zero, corresponding to subsystems that shall be traced over.
1966 The value :math:`-1` refers to the last subsystem.
1967 :type subsystems: int or tuple or list
1969 :param dimensions: Either an integer :math:`d` so that the subsystems
1970 are assumed to be all of shape :math:`d \times d`, or a sequence of
1971 subsystem shapes where an integer :math:`d` within the sequence is
1972 read as :math:`d \times d`. In any case, the elementwise product
1973 over all subsystem shapes must equal the expression's shape.
1974 :type dimensions: int or tuple or list
1976 :raises TypeError: If the subsystems do not match the expression or if
1977 a non-square subsystem is to be traced over.
1978 :raises IndexError: If the subsystem selection is invalid in any other
1979 way.
1981 :Example:
1983 >>> from picos import Constant
1984 >>> A = Constant("A", range(16), (4, 4))
1985 >>> print(A) #doctest: +NORMALIZE_WHITESPACE
1986 [ 0.00e+00 4.00e+00 8.00e+00 1.20e+01]
1987 [ 1.00e+00 5.00e+00 9.00e+00 1.30e+01]
1988 [ 2.00e+00 6.00e+00 1.00e+01 1.40e+01]
1989 [ 3.00e+00 7.00e+00 1.10e+01 1.50e+01]
1990 >>> A0 = A.partial_trace(0); A0
1991 <2×2 Real Constant: A.{tr([2×2])⊗[2×2]}>
1992 >>> print(A0) #doctest: +NORMALIZE_WHITESPACE
1993 [ 1.00e+01 1.80e+01]
1994 [ 1.20e+01 2.00e+01]
1995 >>> A1 = A.partial_trace(1); A1
1996 <2×2 Real Constant: A.{[2×2]⊗tr([2×2])}>
1997 >>> print(A1) #doctest: +NORMALIZE_WHITESPACE
1998 [ 5.00e+00 2.10e+01]
1999 [ 9.00e+00 2.50e+01]
2000 """
2001 # Shape of the original matrix.
2002 m, n = self._shape
2004 if isinstance(dimensions, int):
2005 dimensions = self._square_equal_subsystem_dims(dimensions)
2006 else:
2007 dimensions = [
2008 (d, d) if isinstance(d, int) else d for d in dimensions]
2010 if reduce(
2011 lambda x, y: (x[0]*y[0], x[1]*y[1]), dimensions) != self._shape:
2012 raise TypeError("Subsystem dimensions do not match expression.")
2014 if isinstance(subsystems, int):
2015 subsystems = (subsystems,)
2016 elif not subsystems:
2017 return self
2019 numSys = len(dimensions)
2020 subsystems = set(numSys - 1 if sys == -1 else sys for sys in subsystems)
2022 for sys in subsystems:
2023 if not isinstance(sys, int):
2024 raise IndexError("Subsystem indices must be integer, not {}."
2025 .format(type(sys).__name__))
2026 elif sys < 0:
2027 raise IndexError("Subsystem indices must be nonnegative.")
2028 elif sys >= numSys:
2029 raise IndexError(
2030 "Subsystem index {} out of range for {} systems total."
2031 .format(sys, numSys))
2032 elif dimensions[sys][0] != dimensions[sys][1]:
2033 raise TypeError(
2034 "Subsystem index {} refers to a non-square subsystem that "
2035 "cannot be traced over.".format(sys))
2037 # If all subsystems are traced over, this is the regular trace.
2038 if len(subsystems) == numSys:
2039 return self.tr
2041 # Prepare sparse T such that T·vec(A) = vec(partial_trace(A)).
2042 T = []
2044 # Compute factors of T, one for each subsystem being traced over.
2045 ibh, ibw = m, n
2046 sysStrings = None
2047 for sys in range(numSys):
2048 # Shape of current system.
2049 p, q = dimensions[sys]
2050 sysString = glyphs.matrix(glyphs.shape((p, q)))
2052 # Heigh/width of "outer" blocks whose relative position is
2053 # maintained but that are reduced independently to the size of an
2054 # inner block if the current system is to be traced over.
2055 obh, obw = ibh, ibw
2057 # Height/width of "inner" blocks that are summed if they are
2058 # main-diagonal blocks of an outer block.
2059 ibh, ibw = obh // p, obw // q
2061 # Only trace over selected subsystems.
2062 if sys not in subsystems:
2063 sysStrings = glyphs.kron(sysStrings, sysString) \
2064 if sysStrings else sysString
2065 continue
2066 else:
2067 sysStrings = glyphs.kron(sysStrings, glyphs.trace(sysString)) \
2068 if sysStrings else glyphs.trace(sysString)
2070 # Shape of new matrix.
2071 assert p == q
2072 mN, nN = m // p, n // p
2074 # Prepare one factor of T.
2075 oldLen = m * n
2076 newLen = mN * nN
2077 V, I, J = [1]*(newLen*p), [], []
2078 shape = (newLen, oldLen)
2080 # Determine the summands that make up each entry of the new matrix.
2081 for viN in range(newLen):
2082 # A column/row pair that represents a scalar in the new matrix
2083 # and a number of p scalars within different on-diagonal inner
2084 # blocks but within the same outer block in the old matrix.
2085 cN, rN = divmod(viN, mN)
2087 # Index pair (obi, obj) for outer block in question, row/column
2088 # pair (ibr, ibc) identifies the scalar within each inner block.
2089 obi, ibr = divmod(rN, ibh)
2090 obj, ibc = divmod(cN, ibw)
2092 # Collect summands for the current entry of the new matrix; one
2093 # scalar per on-diagonal inner block.
2094 for k in range(p):
2095 rO = obi*obh + k*ibh + ibr
2096 cO = obj*obw + k*ibw + ibc
2097 I.append(viN)
2098 J.append(cO*m + rO)
2100 # Store the operator that performs the current subsystem trace.
2101 T.insert(0, cvxopt.spmatrix(V, I, J, shape, self._typecode))
2103 # Make next iteration work on the new matrix.
2104 m, n = mN, nN
2106 # Multiply out the linear partial trace operator T.
2107 # TODO: Fast matrix multiplication dynamic program useful here?
2108 T = reduce(lambda A, B: A*B, T)
2110 string = glyphs.ptrace_(self.string, sysStrings)
2111 shape = (m, n)
2112 coefs = {mtbs: T * coef for mtbs, coef in self._coefs.items()}
2114 return self._basetype(string, shape, coefs)
2116 @cached_property
2117 def tr0(self):
2118 r"""Expression with the first :math:`2 \times 2` subsystem traced out.
2120 Only available for a :math:`2^k \times 2^k` matrix with all subsystems
2121 of shape :math:`2 \times 2`. Use :meth:`partial_trace` otherwise.
2122 """
2123 return self.partial_trace(subsystems=0)
2125 @cached_property
2126 def tr1(self):
2127 r"""Expression with the second :math:`2 \times 2` subsystem traced out.
2129 Only available for a :math:`2^k \times 2^k` matrix with all subsystems
2130 of shape :math:`2 \times 2`. Use :meth:`partial_trace` otherwise.
2131 """
2132 return self.partial_trace(subsystems=1)
2134 @cached_property
2135 def tr2(self):
2136 r"""Expression with the third :math:`2 \times 2` subsystem traced out.
2138 Only available for a :math:`2^k \times 2^k` matrix with all subsystems
2139 of shape :math:`2 \times 2`. Use :meth:`partial_trace` otherwise.
2140 """
2141 return self.partial_trace(subsystems=2)
2143 @cached_property
2144 def tr3(self):
2145 r"""Expression with the fourth :math:`2 \times 2` subsystem traced out.
2147 Only available for a :math:`2^k \times 2^k` matrix with all subsystems
2148 of shape :math:`2 \times 2`. Use :meth:`partial_trace` otherwise.
2149 """
2150 return self.partial_trace(subsystems=3)
2152 @cached_property
2153 def trl(self):
2154 r"""Expression with the last :math:`2 \times 2` subsystem traced out.
2156 Only available for a :math:`2^k \times 2^k` matrix with all subsystems
2157 of shape :math:`2 \times 2`. Use :meth:`partial_trace` otherwise.
2158 """
2159 return self.partial_trace(subsystems=-1)
2161 @cached_property
2162 def vec(self):
2163 """Column-major vectorization of the expression as a column vector.
2165 .. note::
2166 Given an expression ``A``, ``A.vec`` and ``A[:]`` produce the same
2167 result (up to its string description) but ``A.vec`` is faster and
2168 its result is cached.
2170 :Example:
2172 >>> from picos import Constant
2173 >>> A = Constant("A", [[1, 2], [3, 4]])
2174 >>> A.vec.equals(A[:])
2175 True
2176 >>> A[:] is A[:]
2177 False
2178 >>> A.vec is A.vec
2179 True
2180 """
2181 if self._shape[1] == 1:
2182 return self
2183 else:
2184 return self._basetype(
2185 glyphs.vec(self.string), (len(self), 1), self._coefs)
2187 def dupvec(self, n):
2188 """Return a (repeated) column-major vectorization of the expression.
2190 :param int n: Number of times to duplicate the vectorization.
2192 :returns: A column vector.
2194 :Example:
2196 >>> from picos import Constant
2197 >>> A = Constant("A", [[1, 2], [3, 4]])
2198 >>> A.dupvec(1) is A.vec
2199 True
2200 >>> A.dupvec(3).equals(A.vec // A.vec // A.vec)
2201 True
2202 """
2203 if not isinstance(n, int):
2204 raise TypeError("Number of copies must be integer.")
2206 if n < 1:
2207 raise ValueError("Number of copies must be positive.")
2209 if n == 1:
2210 return self.vec
2211 else:
2212 string = glyphs.vec(glyphs.comma(self.string, n))
2213 shape = (len(self)*n, 1)
2214 coefs = {mtbs: cvxopt_vcat([coef]*n)
2215 for mtbs, coef in self._coefs.items()}
2217 return self._basetype(string, shape, coefs)
2219 @cached_property
2220 def trilvec(self):
2221 r"""Column-major vectorization of the lower triangular part.
2223 :returns:
2224 A column vector of all elements :math:`A_{ij}` that satisfy
2225 :math:`i \geq j`.
2227 .. note::
2229 If you want a row-major vectorization instead, write ``A.T.triuvec``
2230 instead of ``A.trilvec``.
2232 :Example:
2234 >>> from picos import Constant
2235 >>> A = Constant("A", [[1, 2], [3, 4], [5, 6]])
2236 >>> print(A)
2237 [ 1.00e+00 2.00e+00]
2238 [ 3.00e+00 4.00e+00]
2239 [ 5.00e+00 6.00e+00]
2240 >>> print(A.trilvec)
2241 [ 1.00e+00]
2242 [ 3.00e+00]
2243 [ 5.00e+00]
2244 [ 4.00e+00]
2245 [ 6.00e+00]
2246 """
2247 m, n = self._shape
2249 if n == 1: # Column vector or scalar.
2250 return self
2251 elif m == 1: # Row vector.
2252 return self[0]
2254 # Build a transformation D such that D·vec(A) = trilvec(A).
2255 rows = [j*m + i for j in range(n) for i in range(m) if i >= j]
2256 d = len(rows)
2257 D = cvxopt.spmatrix([1]*d, range(d), rows, (d, len(self)))
2259 string = glyphs.trilvec(self.string)
2260 shape = (d, 1)
2261 coefs = {mtbs: D*coef for mtbs, coef in self._coefs.items()}
2263 return self._basetype(string, shape, coefs)
2265 @cached_property
2266 def triuvec(self):
2267 r"""Column-major vectorization of the upper triangular part.
2269 :returns:
2270 A column vector of all elements :math:`A_{ij}` that satisfy
2271 :math:`i \leq j`.
2273 .. note::
2275 If you want a row-major vectorization instead, write ``A.T.trilvec``
2276 instead of ``A.triuvec``.
2278 :Example:
2280 >>> from picos import Constant
2281 >>> A = Constant("A", [[1, 2, 3], [4, 5, 6]])
2282 >>> print(A)
2283 [ 1.00e+00 2.00e+00 3.00e+00]
2284 [ 4.00e+00 5.00e+00 6.00e+00]
2285 >>> print(A.triuvec)
2286 [ 1.00e+00]
2287 [ 2.00e+00]
2288 [ 5.00e+00]
2289 [ 3.00e+00]
2290 [ 6.00e+00]
2291 """
2292 m, n = self._shape
2294 if m == 1: # Row vector or scalar.
2295 return self
2296 elif n == 1: # Column vector.
2297 return self[0]
2299 # Build a transformation D such that D·vec(A) = triuvec(A).
2300 rows = [j*m + i for j in range(n) for i in range(m) if i <= j]
2301 d = len(rows)
2302 D = cvxopt.spmatrix([1]*d, range(d), rows, (d, len(self)))
2304 string = glyphs.triuvec(self.string)
2305 shape = (d, 1)
2306 coefs = {mtbs: D*coef for mtbs, coef in self._coefs.items()}
2308 return self._basetype(string, shape, coefs)
2310 @cached_property
2311 def svec(self):
2312 """An isometric vectorization of a symmetric or hermitian expression.
2314 In the real symmetric case
2316 - the vectorization format is precisely the one define in [svec]_,
2317 - the vectorization is isometric and isomorphic, and
2318 - this is the same vectorization as used internally by the
2319 :class:`~.variables.SymmetricVariable` class.
2321 In the complex hermitian case
2323 - the same format is used, now resulting in a complex vector,
2324 - the vectorization is isometric but **not** isomorphic as there are
2325 guaranteed zeros in the imaginary part of the vector, and
2326 - this is **not** the same vectorization as the isomorphic,
2327 real-valued one used by :class:`~.variables.HermitianVariable`.
2329 The reverse operation is denoted by :attr:`desvec` in either case.
2331 :raises TypeError:
2332 If the expression is not square.
2334 :raises ValueError:
2335 If the expression is not hermitian.
2336 """
2337 if not self.square:
2338 raise TypeError("Can only compute svec for a square matrix, not for"
2339 " the {} expression {}.".format(self._shape, self.string))
2340 elif not self.hermitian:
2341 raise ValueError("Cannot compute svec for the non-hermitian "
2342 "expression {}.".format(self.string))
2344 vec = SymmetricVectorization(self._shape)
2345 V = vec._full2special
2347 string = glyphs.svec(self.string)
2349 if self.isreal:
2350 vec = SymmetricVectorization(self._shape)
2351 V = vec._full2special
2352 coefs = {mtbs: V*coef for mtbs, coef in self._coefs.items()}
2353 result = self._basetype(string, (vec.dim, 1), coefs)
2354 else:
2355 # NOTE: We need to reproduce svec using triuvec because, for numeric
2356 # reasons, SymmetricVectorization averages the elements in the
2357 # lower and upper triangular part. For symmetric matrices,
2358 # this is equivalent to the formal definition of svec found in
2359 # Datorro's book. For hermitian matrices however it is not.
2360 real_svec = self.real.svec
2361 imag_svec = 2**0.5 * 1j * self.imag.triuvec
2363 result = (real_svec + imag_svec).renamed(string)
2365 with unlocked_cached_properties(result):
2366 result.desvec = self
2368 return result
2370 @cached_property
2371 def desvec(self):
2372 r"""The reverse operation of :attr:`svec`.
2374 :raises TypeError:
2375 If the expression is not a vector or has a dimension outside the
2376 permissible set
2378 .. math::
2380 \left\{ \frac{n(n + 1)}{2}
2381 \mid n \in \mathbb{Z}_{\geq 1} \right\}
2382 = \left\{ n \in \mathbb{Z}_{\geq 1}
2383 \mid \frac{1}{2} \left( \sqrt{8n + 1} - 1 \right)
2384 \in \mathbb{Z}_{\geq 1} \right\}.
2386 :raises ValueError:
2387 In the case of a complex vector, If the vector is not in the image
2388 of :attr:`svec`.
2389 """
2390 if 1 not in self.shape:
2391 raise TypeError(
2392 "Can only compute desvec for a vector, not for the {} "
2393 "expression {}.".format(glyphs.shape(self._shape), self.string))
2395 n = 0.5*((8*len(self) + 1)**0.5 - 1)
2396 if int(n) != n:
2397 raise TypeError("Cannot compute desvec for the {}-dimensional "
2398 "vector {} as this size is not a possible outcome of svec."
2399 .format(len(self), self.string))
2400 n = int(n)
2402 vec = SymmetricVectorization((n, n))
2403 D = vec._special2full
2404 string = glyphs.desvec(self.string)
2406 if self.isreal:
2407 coefs = {mtbs: D*coef for mtbs, coef in self._coefs.items()}
2408 result = self._basetype(string, (n, n), coefs)
2410 assert result.hermitian
2411 else:
2412 # NOTE: While :attr:`svec` performs essentially the same operation
2413 # for both symmetric and hermitian matrices, we now need to
2414 # treat the imaginary separately as the imaginary part of a
2415 # hermitian matrix is skew-symmetric instead of symmetric.
2416 V, I, J = [], D.I, D.J
2417 for v, i, j in zip(D.V, I, J):
2418 V.append(1j*v if i % n <= i // n else -1j*v)
2419 C = cvxopt.spmatrix(V, I, J, D.size, tc="z")
2421 real_desvec = self.real.desvec
2422 imag_desvec = self._basetype(string, (n, n), {
2423 mtbs: C*coef for mtbs, coef in self.imag._coefs.items()})
2425 result = (real_desvec + imag_desvec).renamed(string)
2427 if not result.hermitian:
2428 raise ValueError("The complex vector {} is not in the image of "
2429 "svec. Note that svec is not bijective in the complex case."
2430 .format(self.string))
2432 with unlocked_cached_properties(result):
2433 result.svec = self
2435 return result
2437 def dupdiag(self, n):
2438 """Return a matrix with the (repeated) expression on the diagonal.
2440 Vectorization is performed in column-major order.
2442 :param int n: Number of times to duplicate the vectorization.
2443 """
2444 from .algebra import I
2446 if self.scalar:
2447 return self * I(n)
2449 # Vectorize and duplicate the expression.
2450 vec = self.dupvec(n)
2451 d = len(vec)
2453 # Build a transformation D such that D·vec(A) = vec(diag(vec(A))).
2454 ones = [1]*d
2455 D = cvxopt.spdiag(ones)[:]
2456 D = cvxopt.spmatrix(ones, D.I, range(d), (D.size[0], d))
2458 string = glyphs.diag(vec.string)
2459 shape = (d, d)
2460 coefs = {mtbs: D*coef for mtbs, coef in vec._coefs.items()}
2462 return self._basetype(string, shape, coefs)
2464 @cached_property
2465 def diag(self):
2466 """Diagonal matrix with the expression on the main diagonal.
2468 Vectorization is performed in column-major order.
2469 """
2470 from .algebra import O, I
2472 if self.is0:
2473 return O(len(self), len(self))
2474 elif self.is1:
2475 return I(len(self))
2476 else:
2477 return self.dupdiag(1)
2479 @cached_property
2480 def maindiag(self):
2481 """The main diagonal of the expression as a column vector."""
2482 if 1 in self._shape:
2483 return self[0]
2485 # Build a transformation D such that D·vec(A) = diag(A).
2486 step = self._shape[0] + 1
2487 rows = [i*step for i in range(min(self._shape))]
2488 d = len(rows)
2489 D = cvxopt.spmatrix([1]*d, range(d), rows, (d, len(self)))
2491 string = glyphs.maindiag(self.string)
2492 shape = (d, 1)
2493 coefs = {mtbs: D*coef for mtbs, coef in self._coefs.items()}
2495 return self._basetype(string, shape, coefs)
2497 def factor_out(self, mutable):
2498 r"""Factor out a single mutable from a vector expression.
2500 If this expression is a column vector :math:`a` that depends on some
2501 mutable :math:`x` with a trivial internal vectorization format (i.e.
2502 :class:`~.vectorizations.FullVectorization`), then this method, called
2503 with :math:`x` as its argument, returns a pair of expressions
2504 :math:`(a_x, a_0)` such that :math:`a = a_x\operatorname{vec}(x) + a_0`.
2506 :returns:
2507 Two refined :class:`BiaffineExpression` instances that do not depend
2508 on ``mutable``.
2510 :raises TypeError:
2511 If the expression is not a column vector or if ``mutable`` is not a
2512 :class:`~.mutable.Mutable` or does not have a trivial vectorization
2513 format.
2515 :raises LookupError:
2516 If the expression does not depend on ``mutable``.
2518 :Example:
2520 >>> from picos import RealVariable
2521 >>> from picos.uncertain import UnitBallPerturbationSet
2522 >>> x = RealVariable("x", 3)
2523 >>> z = UnitBallPerturbationSet("z", 3).parameter
2524 >>> a = ((2*x + 3)^z) + 4*x + 5; a
2525 <3×1 Uncertain Affine Expression: (2·x + [3])⊙z + 4·x + [5]>
2526 >>> sorted(a.mutables, key=lambda mtb: mtb.name)
2527 [<3×1 Real Variable: x>, <3×1 Perturbation: z>]
2528 >>> az, a0 = a.factor_out(z)
2529 >>> az
2530 <3×3 Real Affine Expression: ((2·x + [3])⊙z + 4·x + [5])_z>
2531 >>> a0
2532 <3×1 Real Affine Expression: ((2·x + [3])⊙z + 4·x + [5])_0>
2533 >>> sorted(az.mutables.union(a0.mutables), key=lambda mtb: mtb.name)
2534 [<3×1 Real Variable: x>]
2535 >>> (az*z + a0).equals(a)
2536 True
2537 """
2538 if self._shape[1] != 1:
2539 raise TypeError(
2540 "Can only factor out mutables from column vectors but {} has "
2541 "a shape of {}.".format(self.string, glyphs.shape(self._shape)))
2543 mtb = mutable
2545 if not isinstance(mtb, Mutable):
2546 raise TypeError("Can only factor out mutables, not instances of {}."
2547 .format(type(mtb).__name__))
2549 if not isinstance(mtb._vec, FullVectorization):
2550 raise TypeError("Can only factor out mutables with a trivial "
2551 "vectorization format but {} uses {}."
2552 .format(mtb.name, type(mtb._vec).__name__))
2554 if mtb not in self.mutables:
2555 raise LookupError("Cannot factor out {} from {} as the latter does "
2556 "not depend on the former.".format(mtb.name, self.string))
2558 ax_string = glyphs.index(self.string, mtb.name)
2559 ax_shape = (self._shape[0], mtb.dim)
2561 # Recipe in "Robust conic optimization in Python" (Stahlberg 2020).
2562 ax_coefs = {}
2563 for mtbs, coef in self._coefs.items():
2564 if mtb not in set(mtbs):
2565 continue
2566 elif len(mtbs) == 1:
2567 assert mtbs[0] is mtb
2568 self._update_coefs(ax_coefs, (), coef[:])
2569 else:
2570 if mtbs[0] is mtb:
2571 other_mtb = mtbs[1]
2572 C = coef*cvxopt_K(other_mtb.dim, mtb.dim)
2573 else:
2574 assert mtbs[1] is mtb
2575 other_mtb = mtbs[0]
2576 C = coef
2578 d = other_mtb.dim
2579 C = cvxopt.sparse([C[:, i*d:(i+1)*d] for i in range(mtb.dim)])
2580 self._update_coefs(ax_coefs, (other_mtb,), C)
2582 ax = self._basetype(ax_string, ax_shape, ax_coefs)
2584 a0_string = glyphs.index(self.string, 0)
2585 a0_shape = self._shape
2586 a0_coefs = {M: C for M, C in self._coefs.items() if mtb not in set(M)}
2587 a0 = self._basetype(a0_string, a0_shape, a0_coefs)
2589 return ax.refined, a0.refined
2591 # --------------------------------------------------------------------------
2592 # Backwards compatibility methods.
2593 # --------------------------------------------------------------------------
2595 @classmethod
2596 @deprecated("2.0", useInstead="from_constant")
2597 def fromScalar(cls, scalar):
2598 """Create a class instance from a numeric scalar."""
2599 return cls.from_constant(scalar, (1, 1))
2601 @classmethod
2602 @deprecated("2.0", useInstead="from_constant")
2603 def fromMatrix(cls, matrix, size=None):
2604 """Create a class instance from a numeric matrix."""
2605 return cls.from_constant(matrix, size)
2607 @deprecated("2.0", useInstead="object.__xor__")
2608 def hadamard(self, fact):
2609 """Denote the elementwise (or Hadamard) product."""
2610 return self.__xor__(fact)
2612 @deprecated("2.0", useInstead="~.expression.Expression.constant")
2613 def isconstant(self):
2614 """Whether the expression involves no mutables."""
2615 return self.constant
2617 @deprecated("2.0", useInstead="equals")
2618 def same_as(self, other):
2619 """Check mathematical equality with another affine expression."""
2620 return self.equals(other)
2622 @deprecated("2.0", useInstead="T")
2623 def transpose(self):
2624 """Return the matrix transpose."""
2625 return self.T
2627 @cached_property
2628 @deprecated("2.0", useInstead="partial_transpose", decoratorLevel=1)
2629 def Tx(self):
2630 """Auto-detect few subsystems of same shape and transpose the last."""
2631 m, n = self._shape
2632 dims = None
2634 for k in range(2, int(math.log(min(m, n), 2)) + 1):
2635 p, q = int(round(m**(1.0/k))), int(round(n**(1.0/k)))
2636 if m == p**k and n == q**k:
2637 dims = ((p, q),)*k
2638 break
2640 if dims:
2641 return self.partial_transpose(subsystems=-1, dimensions=dims)
2642 else:
2643 raise RuntimeError("Failed to auto-detect subsystem dimensions for "
2644 "partial transposition: Only supported for {} matrices, {}."
2645 .format(glyphs.shape(
2646 (glyphs.power("m", "k"), glyphs.power("n", "k"))),
2647 glyphs.ge("k", 2)))
2649 @deprecated("2.0", useInstead="conj")
2650 def conjugate(self):
2651 """Return the complex conjugate."""
2652 return self.conj
2654 @deprecated("2.0", useInstead="H")
2655 def Htranspose(self):
2656 """Return the conjugate (or Hermitian) transpose."""
2657 return self.H
2659 @deprecated("2.0", reason="PICOS expressions are now immutable.")
2660 def copy(self):
2661 """Return a deep copy of the expression."""
2662 from copy import copy as cp
2664 return self._basetype(cp(self._symbStr), self._shape,
2665 {mtbs: cp(coef) for mtbs, coef in self._coefs.items()})
2667 @deprecated("2.0", reason="PICOS expressions are now immutable.")
2668 def soft_copy(self):
2669 """Return a shallow copy of the expression."""
2670 return self._basetype(self._symbStr, self._shape, self._coefs)
2673# --------------------------------------
2674__all__ = api_end(_API_START, globals())