Coverage for picos/expressions/exp_biaffine.py : 81.42%
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):
1373 """Return the expression reshaped in column-major order.
1375 :Example:
1377 >>> from picos import Constant
1378 >>> C = Constant("C", range(6), (2, 3))
1379 >>> print(C)
1380 [ 0.00e+00 2.00e+00 4.00e+00]
1381 [ 1.00e+00 3.00e+00 5.00e+00]
1382 >>> print(C.reshaped((3, 2)))
1383 [ 0.00e+00 3.00e+00]
1384 [ 1.00e+00 4.00e+00]
1385 [ 2.00e+00 5.00e+00]
1386 """
1387 shape = load_shape(shape, wildcards=True)
1389 if shape == self._shape:
1390 return self
1391 elif shape == (None, None):
1392 return self
1394 length = len(self)
1396 if shape[0] is None:
1397 shape = (length // shape[1], shape[1])
1398 elif shape[1] is None:
1399 shape = (shape[0], length // shape[0])
1401 if shape[0]*shape[1] != length:
1402 raise ValueError("Can only reshape to a matrix of same size.")
1404 string = glyphs.reshaped(self.string, glyphs.shape(shape))
1405 return self._basetype(string, shape, self._coefs)
1407 def broadcasted(self, shape):
1408 """Return the expression broadcasted to the given shape.
1410 :Example:
1412 >>> from picos import Constant
1413 >>> C = Constant("C", range(6), (2, 3))
1414 >>> print(C)
1415 [ 0.00e+00 2.00e+00 4.00e+00]
1416 [ 1.00e+00 3.00e+00 5.00e+00]
1417 >>> print(C.broadcasted((6, 6)))
1418 [ 0.00e+00 2.00e+00 4.00e+00 0.00e+00 2.00e+00 4.00e+00]
1419 [ 1.00e+00 3.00e+00 5.00e+00 1.00e+00 3.00e+00 5.00e+00]
1420 [ 0.00e+00 2.00e+00 4.00e+00 0.00e+00 2.00e+00 4.00e+00]
1421 [ 1.00e+00 3.00e+00 5.00e+00 1.00e+00 3.00e+00 5.00e+00]
1422 [ 0.00e+00 2.00e+00 4.00e+00 0.00e+00 2.00e+00 4.00e+00]
1423 [ 1.00e+00 3.00e+00 5.00e+00 1.00e+00 3.00e+00 5.00e+00]
1424 """
1425 shape = load_shape(shape, wildcards=True)
1426 shape = blend_shapes(shape, self._shape)
1428 if shape == self._shape:
1429 return self
1431 vdup = shape[0] // self._shape[0]
1432 hdup = shape[1] // self._shape[1]
1434 if (self._shape[0] * vdup, self._shape[1] * hdup) != shape:
1435 raise ValueError("Cannot broadcast from shape {} to {}."
1436 .format(glyphs.shape(self._shape), glyphs.shape(shape)))
1438 if self._shape == (1, 1):
1439 string = glyphs.matrix(self.string)
1440 return (self * cvxopt.matrix(1.0, shape)).renamed(string)
1442 string = glyphs.bcasted(self.string, glyphs.shape(shape))
1443 return (cvxopt.matrix(1.0, (vdup, hdup)) @ self).renamed(string)
1445 def reshaped_or_broadcasted(self, shape):
1446 """Return the expression :meth:`reshaped` or :meth:`broadcasted`.
1448 Unlike with :meth:`reshaped` and :meth:`broadcasted`, the target shape
1449 may not contain a wildcard character.
1451 If the wildcard-free target shape has the same number of elements as
1452 the current shape, then this is the same as :meth:`reshaped`, otherwise
1453 it is the same as :meth:`broadcasted`.
1454 """
1455 shape = load_shape(shape, wildcards=False)
1457 try:
1458 if shape[0]*shape[1] == len(self):
1459 return self.reshaped(shape)
1460 else:
1461 return self.broadcasted(shape)
1462 except ValueError:
1463 raise ValueError(
1464 "Cannot reshape or broadcast from shape {} to {}.".format(
1465 glyphs.shape(self._shape), glyphs.shape(shape))) from None
1467 @cached_property
1468 def hermitianized(self):
1469 r"""The expression projected onto the subspace of hermitian matrices.
1471 For a square (complex) affine expression :math:`A`, this is
1472 :math:`\frac{1}{2}(A + A^H)`.
1474 If the expression is not complex, then this is a projection onto the
1475 subspace of symmetric matrices.
1476 """
1477 if not self.square:
1478 raise TypeError("Cannot hermitianize non-square {}.".format(self))
1480 return (self + self.H)/2
1482 @cached_property
1483 def real(self):
1484 """Real part of the expression."""
1485 return self._basetype(glyphs.real(self.string), self._shape,
1486 {mtbs: coef.real() for mtbs, coef in self._coefs.items()})
1488 @cached_property
1489 def imag(self):
1490 """Imaginary part of the expression."""
1491 return self._basetype(glyphs.imag(self.string), self._shape,
1492 {mtbs: coef.imag() for mtbs, coef in self._coefs.items()})
1494 @cached_property
1495 def bilin(self):
1496 """Pure bilinear part of the expression."""
1497 return self._basetype(
1498 glyphs.blinpart(self._symbStr), self._shape, self._bilinear_coefs)
1500 @cached_property
1501 def lin(self):
1502 """Linear part of the expression."""
1503 return self._basetype(
1504 glyphs.linpart(self._symbStr), self._shape, self._linear_coefs)
1506 @cached_property
1507 def noncst(self):
1508 """Nonconstant part of the expression."""
1509 coefs = {mtbs: coefs for mtbs, coefs in self._coefs.items() if mtbs}
1511 return self._basetype(
1512 glyphs.ncstpart(self._symbStr), self._shape, coefs)
1514 @cached_property
1515 def cst(self):
1516 """Constant part of the expression."""
1517 coefs = {(): self._coefs[()]} if () in self._coefs else {}
1519 return self._basetype(glyphs.cstpart(self._symbStr), self._shape, coefs)
1521 @cached_selfinverse_property
1522 def T(self):
1523 """Matrix transpose."""
1524 if len(self) == 1:
1525 return self
1527 m, n = self._shape
1528 K = cvxopt_K(m, n, self._typecode)
1530 string = glyphs.transp(self.string)
1531 shape = (self._shape[1], self._shape[0])
1532 coefs = {mtbs: K * coef for mtbs, coef in self._coefs.items()}
1534 return self._basetype(string, shape, coefs)
1536 @cached_selfinverse_property
1537 def conj(self):
1538 """Complex conjugate."""
1539 string = glyphs.conj(self.string)
1540 coefs = {mtbs: coef.H.T for mtbs, coef in self._coefs.items()}
1542 return self._basetype(string, self._shape, coefs)
1544 @cached_selfinverse_property
1545 def H(self):
1546 """Conjugate (or Hermitian) transpose."""
1547 return self.T.conj.renamed(glyphs.htransp(self._symbStr))
1549 def _square_equal_subsystem_dims(self, diagLen):
1550 """Support :func:`partial_trace` and :func:`partial_transpose`."""
1551 m, n = self._shape
1552 k = math.log(m, diagLen)
1554 if m != n or int(k) != k:
1555 raise TypeError("The expression has shape {} so it cannot be "
1556 "decomposed into subsystems of shape {}.".format(
1557 glyphs.shape(self._shape), glyphs.shape((diagLen,)*2)))
1559 return ((diagLen,)*2,)*int(k)
1561 def partial_transpose(self, subsystems, dimensions=2):
1562 r"""Return the expression with selected subsystems transposed.
1564 If the expression can be written as
1565 :math:`A_0 \otimes \cdots \otimes A_{n-1}` for matrices
1566 :math:`A_0, \ldots, A_{n-1}` with shapes given in ``dimensions``, then
1567 this returns :math:`B_0 \otimes \cdots \otimes B_{n-1}` with
1568 :math:`B_i = A_i^T`, if ``i in subsystems`` (with :math:`i = -1` read as
1569 :math:`n-1`), and :math:`B_i = A_i`, otherwise.
1571 :param subsystems: A collection of or a single subystem number, indexed
1572 from zero, corresponding to subsystems that shall be transposed.
1573 The value :math:`-1` refers to the last subsystem.
1574 :type subsystems: int or tuple or list
1576 :param dimensions: Either an integer :math:`d` so that the subsystems
1577 are assumed to be all of shape :math:`d \times d`, or a sequence of
1578 subsystem shapes where an integer :math:`d` within the sequence is
1579 read as :math:`d \times d`. In any case, the elementwise product
1580 over all subsystem shapes must equal the expression's shape.
1581 :type dimensions: int or tuple or list
1583 :raises TypeError: If the subsystems do not match the expression.
1584 :raises IndexError: If the subsystem selection is invalid.
1586 :Example:
1588 >>> from picos import Constant
1589 >>> A = Constant("A", range(16), (4, 4))
1590 >>> print(A) #doctest: +NORMALIZE_WHITESPACE
1591 [ 0.00e+00 4.00e+00 8.00e+00 1.20e+01]
1592 [ 1.00e+00 5.00e+00 9.00e+00 1.30e+01]
1593 [ 2.00e+00 6.00e+00 1.00e+01 1.40e+01]
1594 [ 3.00e+00 7.00e+00 1.10e+01 1.50e+01]
1595 >>> A0 = A.partial_transpose(0); A0
1596 <4×4 Real Constant: A.{[2×2]ᵀ⊗[2×2]}>
1597 >>> print(A0) #doctest: +NORMALIZE_WHITESPACE
1598 [ 0.00e+00 4.00e+00 2.00e+00 6.00e+00]
1599 [ 1.00e+00 5.00e+00 3.00e+00 7.00e+00]
1600 [ 8.00e+00 1.20e+01 1.00e+01 1.40e+01]
1601 [ 9.00e+00 1.30e+01 1.10e+01 1.50e+01]
1602 >>> A1 = A.partial_transpose(1); A1
1603 <4×4 Real Constant: A.{[2×2]⊗[2×2]ᵀ}>
1604 >>> print(A1) #doctest: +NORMALIZE_WHITESPACE
1605 [ 0.00e+00 1.00e+00 8.00e+00 9.00e+00]
1606 [ 4.00e+00 5.00e+00 1.20e+01 1.30e+01]
1607 [ 2.00e+00 3.00e+00 1.00e+01 1.10e+01]
1608 [ 6.00e+00 7.00e+00 1.40e+01 1.50e+01]
1609 """
1610 m, n = self._shape
1612 if isinstance(dimensions, int):
1613 dimensions = self._square_equal_subsystem_dims(dimensions)
1614 else:
1615 dimensions = [
1616 (d, d) if isinstance(d, int) else d for d in dimensions]
1618 if reduce(
1619 lambda x, y: (x[0]*y[0], x[1]*y[1]), dimensions) != self._shape:
1620 raise TypeError("Subsystem dimensions do not match expression.")
1622 if isinstance(subsystems, int):
1623 subsystems = (subsystems,)
1624 elif not subsystems:
1625 return self
1627 numSys = len(dimensions)
1628 subsystems = set(numSys - 1 if sys == -1 else sys for sys in subsystems)
1630 for sys in subsystems:
1631 if not isinstance(sys, int):
1632 raise IndexError("Subsystem indices must be integer, not {}."
1633 .format(type(sys).__name__))
1634 elif sys < 0:
1635 raise IndexError("Subsystem indices must be nonnegative.")
1636 elif sys >= numSys:
1637 raise IndexError(
1638 "Subsystem index {} out of range for {} systems total."
1639 .format(sys, numSys))
1641 # If all subsystems are transposed, this is regular transposition.
1642 if len(subsystems) == numSys:
1643 return self.T
1645 # Prepare sparse K such that K·vec(A) = vec(partial_transpose(A)).
1646 d = m * n
1647 V = [1]*d
1648 I = range(d)
1649 J = cvxopt.matrix(I)
1650 T = cvxopt.matrix(0, J.size)
1651 obh, obw = 1, 1
1652 sysStrings = None
1654 # Apply transpositions starting with the rightmost Kronecker factor.
1655 for sys in range(numSys - 1, -1, -1):
1656 # Shape of current system.
1657 p, q = dimensions[sys]
1658 sysString = glyphs.matrix(glyphs.shape((p, q)))
1660 # Height/width of "inner" blocks being moved, initially scalars.
1661 ibh, ibw = obh, obw
1663 # Heigh/width of "outer" blocks whose relative position is
1664 # maintained but that are subject to transposition independently.
1665 # In the last iteration this is the shape of the resulting matrix.
1666 obh *= p
1667 obw *= q
1669 # Only transpose selected subsystems.
1670 if sys not in subsystems:
1671 sysStrings = glyphs.kron(sysString, sysStrings) \
1672 if sysStrings else sysString
1673 continue
1674 else:
1675 sysStrings = glyphs.kron(glyphs.transp(sysString), sysStrings) \
1676 if sysStrings else glyphs.transp(sysString)
1678 # Shape of outer blocks after transposition.
1679 obhT, obwT = obw // ibw * ibh, obh // ibh * ibw
1681 # Shape of full matrix after transposition.
1682 mT, nT = m // obh * obhT, n // obw * obwT
1684 for vi in I:
1685 # Full matrix column and row.
1686 c, r = divmod(vi, m)
1688 # Outer block vertical index and row within outer block,
1689 # outer block horizontal index and column within outer block.
1690 obi, obr = divmod(r, obh)
1691 obj, obc = divmod(c, obw)
1693 # Inner block vertical index and row within inner block,
1694 # inner block horizontal index and column within inner block.
1695 ibi, ibr = divmod(obr, ibh)
1696 ibj, ibc = divmod(obc, ibw)
1698 # (1) ibi*ibw + ibc is column within the transposed outer block;
1699 # adding obj*obwT yields the column in the transposed matrix.
1700 # (2) ibj*ibh + ibr is row within the transposed outer block;
1701 # adding obi*obhT yields the row in the transposed matrix.
1702 # (3) tvi is index within the vectorized transposed matrix.
1703 tvi = (obj*obwT + ibi*ibw + ibc)*mT \
1704 + (obi*obhT + ibj*ibh + ibr)
1706 # Prepare the transposition.
1707 T[tvi] = J[vi]
1709 # Apply the transposition.
1710 J, T = T, J
1711 m, n, obh, obw = mT, nT, obhT, obwT
1713 # Finalize the partial commutation matrix.
1714 K = cvxopt.spmatrix(V, I, J, (d, d), self._typecode)
1716 string = glyphs.ptransp_(self.string, sysStrings)
1717 shape = (m, n)
1718 coefs = {mtbs: K * coef for mtbs, coef in self._coefs.items()}
1720 return self._basetype(string, shape, coefs)
1722 @cached_property
1723 def T0(self):
1724 r"""Expression with the first :math:`2 \times 2` subsystem transposed.
1726 Only available for a :math:`2^k \times 2^k` matrix with all subsystems
1727 of shape :math:`2 \times 2`. Use :meth:`partial_transpose` otherwise.
1728 """
1729 return self.partial_transpose(subsystems=0)
1731 @cached_property
1732 def T1(self):
1733 r"""Expression with the second :math:`2 \times 2` subsystem transposed.
1735 Only available for a :math:`2^k \times 2^k` matrix with all subsystems
1736 of shape :math:`2 \times 2`. Use :meth:`partial_transpose` otherwise.
1737 """
1738 return self.partial_transpose(subsystems=1)
1740 @cached_property
1741 def T2(self):
1742 r"""Expression with the third :math:`2 \times 2` subsystem transposed.
1744 Only available for a :math:`2^k \times 2^k` matrix with all subsystems
1745 of shape :math:`2 \times 2`. Use :meth:`partial_transpose` otherwise.
1746 """
1747 return self.partial_transpose(subsystems=2)
1749 @cached_property
1750 def T3(self):
1751 r"""Expression with the fourth :math:`2 \times 2` subsystem transposed.
1753 Only available for a :math:`2^k \times 2^k` matrix with all subsystems
1754 of shape :math:`2 \times 2`. Use :meth:`partial_transpose` otherwise.
1755 """
1756 return self.partial_transpose(subsystems=3)
1758 @cached_property
1759 def Tl(self):
1760 r"""Expression with the last :math:`2 \times 2` subsystem transposed.
1762 Only available for a :math:`2^k \times 2^k` matrix with all subsystems
1763 of shape :math:`2 \times 2`. Use :meth:`partial_transpose` otherwise.
1764 """
1765 return self.partial_transpose(subsystems=-1)
1767 @staticmethod
1768 def _reindex_F(indices, source, destination):
1769 """Convert indices between different tensor shapes in Fortran-order."""
1770 new = []
1771 offset = 0
1772 factor = 1
1774 for index, base in zip(indices, source):
1775 offset += factor*index
1776 factor *= base
1778 for base in destination:
1779 offset, remainder = divmod(offset, base)
1780 new.append(remainder)
1782 return tuple(new)
1784 @staticmethod
1785 def _reindex_C(indices, source, destination):
1786 """Convert indices between different tensor shapes in C-order."""
1787 new = []
1788 offset = 0
1789 factor = 1
1791 for index, base in zip(reversed(indices), reversed(source)):
1792 offset += factor*index
1793 factor *= base
1795 for base in reversed(destination):
1796 offset, remainder = divmod(offset, base)
1797 new.insert(0, remainder)
1799 return tuple(new)
1801 def reshuffled(self, permutation="ikjl", dimensions=None, order="C"):
1802 """Return the reshuffled or realigned expression.
1804 This operation works directly on matrices. However, it is equivalent to
1805 the following sequence of operations:
1807 1. The matrix is reshaped to a tensor with the given ``dimensions`` and
1808 according to ``order``.
1809 2. The tensor's axes are permuted according to ``permutation``.
1810 3. The tensor is reshaped back to the shape of the original matrix
1811 according to ``order``.
1813 For comparison, the following function applies the same operation to a
1814 2D NumPy :class:`~numpy:numpy.ndarray`:
1816 .. code::
1818 def reshuffle_numpy(matrix, permutation, dimensions, order):
1819 P = "{} -> {}".format("".join(sorted(permutation)), permutation)
1820 reshuffled = numpy.reshape(matrix, dimensions, order)
1821 reshuffled = numpy.einsum(P, reshuffled)
1822 return numpy.reshape(reshuffled, matrix.shape, order)
1824 :param permutation:
1825 A sequence of comparable elements with length equal to the number of
1826 tensor dimensions. The sequence is compared to its ordered version
1827 and the resulting permutation pattern is used to permute the tensor
1828 indices. For instance, the string ``"ikjl"`` is compared to its
1829 sorted version ``"ijkl"`` and denotes that the second and third axis
1830 should be swapped.
1831 :type permutation:
1832 str or tuple or list
1834 :param dimensions:
1835 If this is an integer sequence, then it defines the dimensions of
1836 the tensor. If this is :obj:`None`, then the tensor is assumed to be
1837 hypercubic and the number of dimensions is inferred from the
1838 ``permutation`` argument.
1839 :type dimensions:
1840 None or tuple or list
1842 :param str order:
1843 The indexing order to use for the virtual reshaping. Must be either
1844 ``"F"`` for Fortran-order (generalization of column-major) or
1845 ``"C"`` for C-order (generalization of row-major). Note that PICOS
1846 usually reshapes in Fortran-order while NumPy defaults to C-order.
1848 :Example:
1850 >>> from picos import Constant
1851 >>> A = Constant("A", range(16), (4, 4))
1852 >>> print(A) #doctest: +NORMALIZE_WHITESPACE
1853 [ 0.00e+00 4.00e+00 8.00e+00 1.20e+01]
1854 [ 1.00e+00 5.00e+00 9.00e+00 1.30e+01]
1855 [ 2.00e+00 6.00e+00 1.00e+01 1.40e+01]
1856 [ 3.00e+00 7.00e+00 1.10e+01 1.50e+01]
1857 >>> R = A.reshuffled(); R
1858 <4×4 Real Constant: shuffled(A,ikjl,C)>
1859 >>> print(R) #doctest: +NORMALIZE_WHITESPACE
1860 [ 0.00e+00 4.00e+00 1.00e+00 5.00e+00]
1861 [ 8.00e+00 1.20e+01 9.00e+00 1.30e+01]
1862 [ 2.00e+00 6.00e+00 3.00e+00 7.00e+00]
1863 [ 1.00e+01 1.40e+01 1.10e+01 1.50e+01]
1864 >>> A.reshuffled("ji").equals(A.T) # Regular transposition.
1865 True
1866 >>> A.reshuffled("3214").equals(A.T0) # Partial transposition (1).
1867 True
1868 >>> A.reshuffled("1432").equals(A.T1) # Partial transposition (2).
1869 True
1870 """
1871 m, n = self._shape
1872 mn = m*n
1874 # Load the permutation.
1875 ordered = sorted(permutation)
1876 P = dict(enumerate(ordered.index(element) for element in permutation))
1878 if len(set(P.values())) < len(P):
1879 raise ValueError("The sequence defining the permutation appears to "
1880 "contain duplicate elements.")
1882 assert not set(P.keys()).symmetric_difference(set(P.values()))
1884 numDims = len(P)
1886 # Load the dimensions.
1887 guessDimensions = dimensions is None
1889 if guessDimensions:
1890 dimensions = (int(mn**(1.0 / numDims)),)*numDims
1891 else:
1892 if len(dimensions) != numDims:
1893 raise ValueError("The number of indices does not match the "
1894 "number of dimensions.")
1896 if reduce(int.__mul__, dimensions, 1) != mn:
1897 raise TypeError("The {} matrix {} cannot be reshaped to a {} "
1898 "tensor.".format(glyphs.shape(self.shape), self.string,
1899 "hypercubic order {}".format(numDims) if guessDimensions
1900 else glyphs.size("", "").join(str(d) for d in dimensions)))
1902 # Load the indexing order.
1903 if order not in "FC":
1904 raise ValueError("Order must be given as 'F' or 'C'.")
1906 reindex = self._reindex_F if order == "F" else self._reindex_C
1908 # Nothing to do for the neutral permutation.
1909 if all(key == val for key, val in P.items()):
1910 return self
1912 # Create a sparse mtrix R such that R·vec(A) = vec(reshuffle(A)).
1913 V, I, J = [1]*mn, [], range(mn)
1914 for i in range(mn):
1915 (k, j) = divmod(i, m) # (j, k) are column-major matrix indices.
1916 indices = reindex((j, k), (m, n), dimensions)
1917 newIndices = tuple(indices[P[d]] for d in range(numDims))
1918 newDimensions = tuple(dimensions[P[d]] for d in range(numDims))
1919 (j, k) = reindex(newIndices, newDimensions, (m, n))
1920 I.append(k*m + j)
1921 R = cvxopt.spmatrix(V, I, J, (mn, mn), self._typecode)
1923 # Create the string.
1924 strArgs = [self.string, str(permutation).replace(" ", ""), order]
1926 if not guessDimensions:
1927 strArgs.insert(2, str(dimensions).replace(" ", ""))
1929 string = glyphs.shuffled(",".join(strArgs))
1931 # Finalize the new expression.
1932 shape = (m, n)
1933 coefs = {mtbs: R * coef for mtbs, coef in self._coefs.items()}
1935 return self._basetype(string, shape, coefs)
1937 @cached_property
1938 def sum(self):
1939 """Sum over all scalar elements of the expression."""
1940 # NOTE: glyphs.clever_dotp detects this case and uses the sum glyph.
1941 # NOTE: 1 on the right hand side in case self is complex.
1942 return (self | 1)
1944 @cached_property
1945 def tr(self):
1946 """Trace of a square expression."""
1947 if not self.square:
1948 raise TypeError("Cannot compute the trace of non-square {}."
1949 .format(self.string))
1951 # NOTE: glyphs.clever_dotp detects this case and uses the trace glyph.
1952 # NOTE: "I" on the right hand side in case self is complex.
1953 return (self | "I")
1955 def partial_trace(self, subsystems, dimensions=2):
1956 r"""Return the partial trace over selected subsystems.
1958 If the expression can be written as
1959 :math:`A_0 \otimes \cdots \otimes A_{n-1}` for matrices
1960 :math:`A_0, \ldots, A_{n-1}` with shapes given in ``dimensions``, then
1961 this returns :math:`B_0 \otimes \cdots \otimes B_{n-1}` with
1962 :math:`B_i = \operatorname{tr}(A_i)`, if ``i in subsystems`` (with
1963 :math:`i = -1` read as :math:`n-1`), and :math:`B_i = A_i`, otherwise.
1965 :param subsystems: A collection of or a single subystem number, indexed
1966 from zero, corresponding to subsystems that shall be traced over.
1967 The value :math:`-1` refers to the last subsystem.
1968 :type subsystems: int or tuple or list
1970 :param dimensions: Either an integer :math:`d` so that the subsystems
1971 are assumed to be all of shape :math:`d \times d`, or a sequence of
1972 subsystem shapes where an integer :math:`d` within the sequence is
1973 read as :math:`d \times d`. In any case, the elementwise product
1974 over all subsystem shapes must equal the expression's shape.
1975 :type dimensions: int or tuple or list
1977 :raises TypeError: If the subsystems do not match the expression or if
1978 a non-square subsystem is to be traced over.
1979 :raises IndexError: If the subsystem selection is invalid in any other
1980 way.
1982 :Example:
1984 >>> from picos import Constant
1985 >>> A = Constant("A", range(16), (4, 4))
1986 >>> print(A) #doctest: +NORMALIZE_WHITESPACE
1987 [ 0.00e+00 4.00e+00 8.00e+00 1.20e+01]
1988 [ 1.00e+00 5.00e+00 9.00e+00 1.30e+01]
1989 [ 2.00e+00 6.00e+00 1.00e+01 1.40e+01]
1990 [ 3.00e+00 7.00e+00 1.10e+01 1.50e+01]
1991 >>> A0 = A.partial_trace(0); A0
1992 <2×2 Real Constant: A.{tr([2×2])⊗[2×2]}>
1993 >>> print(A0) #doctest: +NORMALIZE_WHITESPACE
1994 [ 1.00e+01 1.80e+01]
1995 [ 1.20e+01 2.00e+01]
1996 >>> A1 = A.partial_trace(1); A1
1997 <2×2 Real Constant: A.{[2×2]⊗tr([2×2])}>
1998 >>> print(A1) #doctest: +NORMALIZE_WHITESPACE
1999 [ 5.00e+00 2.10e+01]
2000 [ 9.00e+00 2.50e+01]
2001 """
2002 # Shape of the original matrix.
2003 m, n = self._shape
2005 if isinstance(dimensions, int):
2006 dimensions = self._square_equal_subsystem_dims(dimensions)
2007 else:
2008 dimensions = [
2009 (d, d) if isinstance(d, int) else d for d in dimensions]
2011 if reduce(
2012 lambda x, y: (x[0]*y[0], x[1]*y[1]), dimensions) != self._shape:
2013 raise TypeError("Subsystem dimensions do not match expression.")
2015 if isinstance(subsystems, int):
2016 subsystems = (subsystems,)
2017 elif not subsystems:
2018 return self
2020 numSys = len(dimensions)
2021 subsystems = set(numSys - 1 if sys == -1 else sys for sys in subsystems)
2023 for sys in subsystems:
2024 if not isinstance(sys, int):
2025 raise IndexError("Subsystem indices must be integer, not {}."
2026 .format(type(sys).__name__))
2027 elif sys < 0:
2028 raise IndexError("Subsystem indices must be nonnegative.")
2029 elif sys >= numSys:
2030 raise IndexError(
2031 "Subsystem index {} out of range for {} systems total."
2032 .format(sys, numSys))
2033 elif dimensions[sys][0] != dimensions[sys][1]:
2034 raise TypeError(
2035 "Subsystem index {} refers to a non-square subsystem that "
2036 "cannot be traced over.".format(sys))
2038 # If all subsystems are traced over, this is the regular trace.
2039 if len(subsystems) == numSys:
2040 return self.tr
2042 # Prepare sparse T such that T·vec(A) = vec(partial_trace(A)).
2043 T = []
2045 # Compute factors of T, one for each subsystem being traced over.
2046 ibh, ibw = m, n
2047 sysStrings = None
2048 for sys in range(numSys):
2049 # Shape of current system.
2050 p, q = dimensions[sys]
2051 sysString = glyphs.matrix(glyphs.shape((p, q)))
2053 # Heigh/width of "outer" blocks whose relative position is
2054 # maintained but that are reduced independently to the size of an
2055 # inner block if the current system is to be traced over.
2056 obh, obw = ibh, ibw
2058 # Height/width of "inner" blocks that are summed if they are
2059 # main-diagonal blocks of an outer block.
2060 ibh, ibw = obh // p, obw // q
2062 # Only trace over selected subsystems.
2063 if sys not in subsystems:
2064 sysStrings = glyphs.kron(sysStrings, sysString) \
2065 if sysStrings else sysString
2066 continue
2067 else:
2068 sysStrings = glyphs.kron(sysStrings, glyphs.trace(sysString)) \
2069 if sysStrings else glyphs.trace(sysString)
2071 # Shape of new matrix.
2072 assert p == q
2073 mN, nN = m // p, n // p
2075 # Prepare one factor of T.
2076 oldLen = m * n
2077 newLen = mN * nN
2078 V, I, J = [1]*(newLen*p), [], []
2079 shape = (newLen, oldLen)
2081 # Determine the summands that make up each entry of the new matrix.
2082 for viN in range(newLen):
2083 # A column/row pair that represents a scalar in the new matrix
2084 # and a number of p scalars within different on-diagonal inner
2085 # blocks but within the same outer block in the old matrix.
2086 cN, rN = divmod(viN, mN)
2088 # Index pair (obi, obj) for outer block in question, row/column
2089 # pair (ibr, ibc) identifies the scalar within each inner block.
2090 obi, ibr = divmod(rN, ibh)
2091 obj, ibc = divmod(cN, ibw)
2093 # Collect summands for the current entry of the new matrix; one
2094 # scalar per on-diagonal inner block.
2095 for k in range(p):
2096 rO = obi*obh + k*ibh + ibr
2097 cO = obj*obw + k*ibw + ibc
2098 I.append(viN)
2099 J.append(cO*m + rO)
2101 # Store the operator that performs the current subsystem trace.
2102 T.insert(0, cvxopt.spmatrix(V, I, J, shape, self._typecode))
2104 # Make next iteration work on the new matrix.
2105 m, n = mN, nN
2107 # Multiply out the linear partial trace operator T.
2108 # TODO: Fast matrix multiplication dynamic program useful here?
2109 T = reduce(lambda A, B: A*B, T)
2111 string = glyphs.ptrace_(self.string, sysStrings)
2112 shape = (m, n)
2113 coefs = {mtbs: T * coef for mtbs, coef in self._coefs.items()}
2115 return self._basetype(string, shape, coefs)
2117 @cached_property
2118 def tr0(self):
2119 r"""Expression with the first :math:`2 \times 2` subsystem traced out.
2121 Only available for a :math:`2^k \times 2^k` matrix with all subsystems
2122 of shape :math:`2 \times 2`. Use :meth:`partial_trace` otherwise.
2123 """
2124 return self.partial_trace(subsystems=0)
2126 @cached_property
2127 def tr1(self):
2128 r"""Expression with the second :math:`2 \times 2` subsystem traced out.
2130 Only available for a :math:`2^k \times 2^k` matrix with all subsystems
2131 of shape :math:`2 \times 2`. Use :meth:`partial_trace` otherwise.
2132 """
2133 return self.partial_trace(subsystems=1)
2135 @cached_property
2136 def tr2(self):
2137 r"""Expression with the third :math:`2 \times 2` subsystem traced out.
2139 Only available for a :math:`2^k \times 2^k` matrix with all subsystems
2140 of shape :math:`2 \times 2`. Use :meth:`partial_trace` otherwise.
2141 """
2142 return self.partial_trace(subsystems=2)
2144 @cached_property
2145 def tr3(self):
2146 r"""Expression with the fourth :math:`2 \times 2` subsystem traced out.
2148 Only available for a :math:`2^k \times 2^k` matrix with all subsystems
2149 of shape :math:`2 \times 2`. Use :meth:`partial_trace` otherwise.
2150 """
2151 return self.partial_trace(subsystems=3)
2153 @cached_property
2154 def trl(self):
2155 r"""Expression with the last :math:`2 \times 2` subsystem traced out.
2157 Only available for a :math:`2^k \times 2^k` matrix with all subsystems
2158 of shape :math:`2 \times 2`. Use :meth:`partial_trace` otherwise.
2159 """
2160 return self.partial_trace(subsystems=-1)
2162 @cached_property
2163 def vec(self):
2164 """Column-major vectorization of the expression as a column vector.
2166 .. note::
2167 Given an expression ``A``, ``A.vec`` and ``A[:]`` produce the same
2168 result (up to its string description) but ``A.vec`` is faster and
2169 its result is cached.
2171 :Example:
2173 >>> from picos import Constant
2174 >>> A = Constant("A", [[1, 2], [3, 4]])
2175 >>> A.vec.equals(A[:])
2176 True
2177 >>> A[:] is A[:]
2178 False
2179 >>> A.vec is A.vec
2180 True
2181 """
2182 if self._shape[1] == 1:
2183 return self
2184 else:
2185 return self._basetype(
2186 glyphs.vec(self.string), (len(self), 1), self._coefs)
2188 def dupvec(self, n):
2189 """Return a (repeated) column-major vectorization of the expression.
2191 :param int n: Number of times to duplicate the vectorization.
2193 :returns: A column vector.
2195 :Example:
2197 >>> from picos import Constant
2198 >>> A = Constant("A", [[1, 2], [3, 4]])
2199 >>> A.dupvec(1) is A.vec
2200 True
2201 >>> A.dupvec(3).equals(A.vec // A.vec // A.vec)
2202 True
2203 """
2204 if not isinstance(n, int):
2205 raise TypeError("Number of copies must be integer.")
2207 if n < 1:
2208 raise ValueError("Number of copies must be positive.")
2210 if n == 1:
2211 return self.vec
2212 else:
2213 string = glyphs.vec(glyphs.comma(self.string, n))
2214 shape = (len(self)*n, 1)
2215 coefs = {mtbs: cvxopt_vcat([coef]*n)
2216 for mtbs, coef in self._coefs.items()}
2218 return self._basetype(string, shape, coefs)
2220 @cached_property
2221 def trilvec(self):
2222 r"""Column-major vectorization of the lower triangular part.
2224 :returns:
2225 A column vector of all elements :math:`A_{ij}` that satisfy
2226 :math:`i \geq j`.
2228 .. note::
2230 If you want a row-major vectorization instead, write ``A.T.triuvec``
2231 instead of ``A.trilvec``.
2233 :Example:
2235 >>> from picos import Constant
2236 >>> A = Constant("A", [[1, 2], [3, 4], [5, 6]])
2237 >>> print(A)
2238 [ 1.00e+00 2.00e+00]
2239 [ 3.00e+00 4.00e+00]
2240 [ 5.00e+00 6.00e+00]
2241 >>> print(A.trilvec)
2242 [ 1.00e+00]
2243 [ 3.00e+00]
2244 [ 5.00e+00]
2245 [ 4.00e+00]
2246 [ 6.00e+00]
2247 """
2248 m, n = self._shape
2250 if n == 1: # Column vector or scalar.
2251 return self
2252 elif m == 1: # Row vector.
2253 return self[0]
2255 # Build a transformation D such that D·vec(A) = trilvec(A).
2256 rows = [j*m + i for j in range(n) for i in range(m) if i >= j]
2257 d = len(rows)
2258 D = cvxopt.spmatrix([1]*d, range(d), rows, (d, len(self)))
2260 string = glyphs.trilvec(self.string)
2261 shape = (d, 1)
2262 coefs = {mtbs: D*coef for mtbs, coef in self._coefs.items()}
2264 return self._basetype(string, shape, coefs)
2266 @cached_property
2267 def triuvec(self):
2268 r"""Column-major vectorization of the upper triangular part.
2270 :returns:
2271 A column vector of all elements :math:`A_{ij}` that satisfy
2272 :math:`i \leq j`.
2274 .. note::
2276 If you want a row-major vectorization instead, write ``A.T.trilvec``
2277 instead of ``A.triuvec``.
2279 :Example:
2281 >>> from picos import Constant
2282 >>> A = Constant("A", [[1, 2, 3], [4, 5, 6]])
2283 >>> print(A)
2284 [ 1.00e+00 2.00e+00 3.00e+00]
2285 [ 4.00e+00 5.00e+00 6.00e+00]
2286 >>> print(A.triuvec)
2287 [ 1.00e+00]
2288 [ 2.00e+00]
2289 [ 5.00e+00]
2290 [ 3.00e+00]
2291 [ 6.00e+00]
2292 """
2293 m, n = self._shape
2295 if m == 1: # Row vector or scalar.
2296 return self
2297 elif n == 1: # Column vector.
2298 return self[0]
2300 # Build a transformation D such that D·vec(A) = triuvec(A).
2301 rows = [j*m + i for j in range(n) for i in range(m) if i <= j]
2302 d = len(rows)
2303 D = cvxopt.spmatrix([1]*d, range(d), rows, (d, len(self)))
2305 string = glyphs.triuvec(self.string)
2306 shape = (d, 1)
2307 coefs = {mtbs: D*coef for mtbs, coef in self._coefs.items()}
2309 return self._basetype(string, shape, coefs)
2311 @cached_property
2312 def svec(self):
2313 """An isometric vectorization of a symmetric or hermitian expression.
2315 In the real symmetric case
2317 - the vectorization format is precisely the one define in [svec]_,
2318 - the vectorization is isometric and isomorphic, and
2319 - this is the same vectorization as used internally by the
2320 :class:`~.variables.SymmetricVariable` class.
2322 In the complex hermitian case
2324 - the same format is used, now resulting in a complex vector,
2325 - the vectorization is isometric but **not** isomorphic as there are
2326 guaranteed zeros in the imaginary part of the vector, and
2327 - this is **not** the same vectorization as the isomorphic,
2328 real-valued one used by :class:`~.variables.HermitianVariable`.
2330 The reverse operation is denoted by :attr:`desvec` in either case.
2332 :raises TypeError:
2333 If the expression is not square.
2335 :raises ValueError:
2336 If the expression is not hermitian.
2337 """
2338 if not self.square:
2339 raise TypeError("Can only compute svec for a square matrix, not for"
2340 " the {} expression {}.".format(self._shape, self.string))
2341 elif not self.hermitian:
2342 raise ValueError("Cannot compute svec for the non-hermitian "
2343 "expression {}.".format(self.string))
2345 vec = SymmetricVectorization(self._shape)
2346 V = vec._full2special
2348 string = glyphs.svec(self.string)
2350 if self.isreal:
2351 vec = SymmetricVectorization(self._shape)
2352 V = vec._full2special
2353 coefs = {mtbs: V*coef for mtbs, coef in self._coefs.items()}
2354 result = self._basetype(string, (vec.dim, 1), coefs)
2355 else:
2356 # NOTE: We need to reproduce svec using triuvec because, for numeric
2357 # reasons, SymmetricVectorization averages the elements in the
2358 # lower and upper triangular part. For symmetric matrices,
2359 # this is equivalent to the formal definition of svec found in
2360 # Datorro's book. For hermitian matrices however it is not.
2361 real_svec = self.real.svec
2362 imag_svec = 2**0.5 * 1j * self.imag.triuvec
2364 result = (real_svec + imag_svec).renamed(string)
2366 with unlocked_cached_properties(result):
2367 result.desvec = self
2369 return result
2371 @cached_property
2372 def desvec(self):
2373 r"""The reverse operation of :attr:`svec`.
2375 :raises TypeError:
2376 If the expression is not a vector or has a dimension outside the
2377 permissible set
2379 .. math::
2381 \left\{ \frac{n(n + 1)}{2}
2382 \mid n \in \mathbb{Z}_{\geq 1} \right\}
2383 = \left\{ n \in \mathbb{Z}_{\geq 1}
2384 \mid \frac{1}{2} \left( \sqrt{8n + 1} - 1 \right)
2385 \in \mathbb{Z}_{\geq 1} \right\}.
2387 :raises ValueError:
2388 In the case of a complex vector, If the vector is not in the image
2389 of :attr:`svec`.
2390 """
2391 if 1 not in self.shape:
2392 raise TypeError(
2393 "Can only compute desvec for a vector, not for the {} "
2394 "expression {}.".format(glyphs.shape(self._shape), self.string))
2396 n = 0.5*((8*len(self) + 1)**0.5 - 1)
2397 if int(n) != n:
2398 raise TypeError("Cannot compute desvec for the {}-dimensional "
2399 "vector {} as this size is not a possible outcome of svec."
2400 .format(len(self), self.string))
2401 n = int(n)
2403 vec = SymmetricVectorization((n, n))
2404 D = vec._special2full
2405 string = glyphs.desvec(self.string)
2407 if self.isreal:
2408 coefs = {mtbs: D*coef for mtbs, coef in self._coefs.items()}
2409 result = self._basetype(string, (n, n), coefs)
2411 assert result.hermitian
2412 else:
2413 # NOTE: While :attr:`svec` performs essentially the same operation
2414 # for both symmetric and hermitian matrices, we now need to
2415 # treat the imaginary separately as the imaginary part of a
2416 # hermitian matrix is skew-symmetric instead of symmetric.
2417 V, I, J = [], D.I, D.J
2418 for v, i, j in zip(D.V, I, J):
2419 V.append(1j*v if i % n <= i // n else -1j*v)
2420 C = cvxopt.spmatrix(V, I, J, D.size, tc="z")
2422 real_desvec = self.real.desvec
2423 imag_desvec = self._basetype(string, (n, n), {
2424 mtbs: C*coef for mtbs, coef in self.imag._coefs.items()})
2426 result = (real_desvec + imag_desvec).renamed(string)
2428 if not result.hermitian:
2429 raise ValueError("The complex vector {} is not in the image of "
2430 "svec. Note that svec is not bijective in the complex case."
2431 .format(self.string))
2433 with unlocked_cached_properties(result):
2434 result.svec = self
2436 return result
2438 def dupdiag(self, n):
2439 """Return a matrix with the (repeated) expression on the diagonal.
2441 Vectorization is performed in column-major order.
2443 :param int n: Number of times to duplicate the vectorization.
2444 """
2445 # Vectorize and duplicate the expression.
2446 vec = self.dupvec(n)
2447 d = len(vec)
2449 # Build a transformation D such that D·vec(A) = vec(diag(vec(A))).
2450 ones = [1]*d
2451 D = cvxopt.spdiag(ones)[:]
2452 D = cvxopt.spmatrix(ones, D.I, range(d), (D.size[0], d))
2454 string = glyphs.diag(vec.string)
2455 shape = (d, d)
2456 coefs = {mtbs: D*coef for mtbs, coef in vec._coefs.items()}
2458 return self._basetype(string, shape, coefs)
2460 @cached_property
2461 def diag(self):
2462 """Diagonal matrix with the expression on the main diagonal.
2464 Vectorization is performed in column-major order.
2465 """
2466 return self.dupdiag(1)
2468 @cached_property
2469 def maindiag(self):
2470 """The main diagonal of the expression as a column vector."""
2471 if 1 in self._shape:
2472 return self[0]
2474 # Build a transformation D such that D·vec(A) = diag(A).
2475 step = self._shape[0] + 1
2476 rows = [i*step for i in range(min(self._shape))]
2477 d = len(rows)
2478 D = cvxopt.spmatrix([1]*d, range(d), rows, (d, len(self)))
2480 string = glyphs.maindiag(self.string)
2481 shape = (d, 1)
2482 coefs = {mtbs: D*coef for mtbs, coef in self._coefs.items()}
2484 return self._basetype(string, shape, coefs)
2486 def factor_out(self, mutable):
2487 r"""Factor out a single mutable from a vector expression.
2489 If this expression is a column vector :math:`a` that depends on some
2490 mutable :math:`x` with a trivial internal vectorization format (i.e.
2491 :class:`~.vectorizations.FullVectorization`), then this method, called
2492 with :math:`x` as its argument, returns a pair of expressions
2493 :math:`(a_x, a_0)` such that :math:`a = a_x\operatorname{vec}(x) + a_0`.
2495 :returns:
2496 Two refined :class:`BiaffineExpression` instances that do not depend
2497 on ``mutable``.
2499 :raises TypeError:
2500 If the expression is not a column vector or if ``mutable`` is not a
2501 :class:`~.mutable.Mutable` or does not have a trivial vectorization
2502 format.
2504 :raises LookupError:
2505 If the expression does not depend on ``mutable``.
2507 :Example:
2509 >>> from picos import RealVariable
2510 >>> from picos.uncertain import UnitBallPerturbationSet
2511 >>> x = RealVariable("x", 3)
2512 >>> z = UnitBallPerturbationSet("z", 3).parameter
2513 >>> a = ((2*x + 3)^z) + 4*x + 5; a
2514 <3×1 Uncertain Affine Expression: (2·x + [3])⊙z + 4·x + [5]>
2515 >>> sorted(a.mutables, key=lambda mtb: mtb.name)
2516 [<3×1 Real Variable: x>, <3×1 Perturbation: z>]
2517 >>> az, a0 = a.factor_out(z)
2518 >>> az
2519 <3×3 Real Affine Expression: ((2·x + [3])⊙z + 4·x + [5])_z>
2520 >>> a0
2521 <3×1 Real Affine Expression: ((2·x + [3])⊙z + 4·x + [5])_0>
2522 >>> sorted(az.mutables.union(a0.mutables), key=lambda mtb: mtb.name)
2523 [<3×1 Real Variable: x>]
2524 >>> (az*z + a0).equals(a)
2525 True
2526 """
2527 if self._shape[1] != 1:
2528 raise TypeError(
2529 "Can only factor out mutables from column vectors but {} has "
2530 "a shape of {}.".format(self.string, glyphs.shape(self._shape)))
2532 mtb = mutable
2534 if not isinstance(mtb, Mutable):
2535 raise TypeError("Can only factor out mutables, not instances of {}."
2536 .format(type(mtb).__name__))
2538 if not isinstance(mtb._vec, FullVectorization):
2539 raise TypeError("Can only factor out mutables with a trivial "
2540 "vectorization format but {} uses {}."
2541 .format(mtb.name, type(mtb._vec).__name__))
2543 if mtb not in self.mutables:
2544 raise LookupError("Cannot factor out {} from {} as the latter does "
2545 "not depend on the former.".format(mtb.name, self.string))
2547 ax_string = glyphs.index(self.string, mtb.name)
2548 ax_shape = (self._shape[0], mtb.dim)
2550 # Recipe in "Robust conic optimization in Python" (Stahlberg 2020).
2551 ax_coefs = {}
2552 for mtbs, coef in self._coefs.items():
2553 if mtb not in set(mtbs):
2554 continue
2555 elif len(mtbs) == 1:
2556 assert mtbs[0] is mtb
2557 self._update_coefs(ax_coefs, (), coef[:])
2558 else:
2559 if mtbs[0] is mtb:
2560 other_mtb = mtbs[1]
2561 C = coef*cvxopt_K(other_mtb.dim, mtb.dim)
2562 else:
2563 assert mtbs[1] is mtb
2564 other_mtb = mtbs[0]
2565 C = coef
2567 d = other_mtb.dim
2568 C = cvxopt.sparse([C[:, i*d:(i+1)*d] for i in range(mtb.dim)])
2569 self._update_coefs(ax_coefs, (other_mtb,), C)
2571 ax = self._basetype(ax_string, ax_shape, ax_coefs)
2573 a0_string = glyphs.index(self.string, 0)
2574 a0_shape = self._shape
2575 a0_coefs = {M: C for M, C in self._coefs.items() if mtb not in set(M)}
2576 a0 = self._basetype(a0_string, a0_shape, a0_coefs)
2578 return ax.refined, a0.refined
2580 # --------------------------------------------------------------------------
2581 # Backwards compatibility methods.
2582 # --------------------------------------------------------------------------
2584 @classmethod
2585 @deprecated("2.0", useInstead="from_constant")
2586 def fromScalar(cls, scalar):
2587 """Create a class instance from a numeric scalar."""
2588 return cls.from_constant(scalar, (1, 1))
2590 @classmethod
2591 @deprecated("2.0", useInstead="from_constant")
2592 def fromMatrix(cls, matrix, size=None):
2593 """Create a class instance from a numeric matrix."""
2594 return cls.from_constant(matrix, size)
2596 @deprecated("2.0", useInstead="object.__xor__")
2597 def hadamard(self, fact):
2598 """Denote the elementwise (or Hadamard) product."""
2599 return self.__xor__(fact)
2601 @deprecated("2.0", useInstead="~.expression.Expression.constant")
2602 def isconstant(self):
2603 """Whether the expression involves no mutables."""
2604 return self.constant
2606 @deprecated("2.0", useInstead="equals")
2607 def same_as(self, other):
2608 """Check mathematical equality with another affine expression."""
2609 return self.equals(other)
2611 @deprecated("2.0", useInstead="T")
2612 def transpose(self):
2613 """Return the matrix transpose."""
2614 return self.T
2616 @cached_property
2617 @deprecated("2.0", useInstead="partial_transpose", decoratorLevel=1)
2618 def Tx(self):
2619 """Auto-detect few subsystems of same shape and transpose the last."""
2620 m, n = self._shape
2621 dims = None
2623 for k in range(2, int(math.log(min(m, n), 2)) + 1):
2624 p, q = int(round(m**(1.0/k))), int(round(n**(1.0/k)))
2625 if m == p**k and n == q**k:
2626 dims = ((p, q),)*k
2627 break
2629 if dims:
2630 return self.partial_transpose(subsystems=-1, dimensions=dims)
2631 else:
2632 raise RuntimeError("Failed to auto-detect subsystem dimensions for "
2633 "partial transposition: Only supported for {} matrices, {}."
2634 .format(glyphs.shape(
2635 (glyphs.power("m", "k"), glyphs.power("n", "k"))),
2636 glyphs.ge("k", 2)))
2638 @deprecated("2.0", useInstead="conj")
2639 def conjugate(self):
2640 """Return the complex conjugate."""
2641 return self.conj
2643 @deprecated("2.0", useInstead="H")
2644 def Htranspose(self):
2645 """Return the conjugate (or Hermitian) transpose."""
2646 return self.H
2648 @deprecated("2.0", reason="PICOS expressions are now immutable.")
2649 def copy(self):
2650 """Return a deep copy of the expression."""
2651 from copy import copy as cp
2653 return self._basetype(cp(self._symbStr), self._shape,
2654 {mtbs: cp(coef) for mtbs, coef in self._coefs.items()})
2656 @deprecated("2.0", reason="PICOS expressions are now immutable.")
2657 def soft_copy(self):
2658 """Return a shallow copy of the expression."""
2659 return self._basetype(self._symbStr, self._shape, self._coefs)
2662# --------------------------------------
2663__all__ = api_end(_API_START, globals())