# coding: utf-8
# ------------------------------------------------------------------------------
# Copyright (C) 2019 Maximilian Stahlberg
# Based on the original picos.expressions module by Guillaume Sagnol.
#
# This file is part of PICOS.
#
# PICOS is free software: you can redistribute it and/or modify
# it under the terms of the GNU General Public License as published by
# the Free Software Foundation, either version 3 of the License, or
# (at your option) any later version.
#
# PICOS is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
# GNU General Public License for more details.
#
# You should have received a copy of the GNU General Public License
# along with this program. If not, see <http://www.gnu.org/licenses/>.
# ------------------------------------------------------------------------------
"""Implements affine expression types."""
import math
import operator
from collections import namedtuple
from functools import reduce
import cvxopt
import numpy
from .. import glyphs
from ..apidoc import api_end, api_start
from ..caching import (cached_property, cached_selfinverse_property,
cached_selfinverse_unary_operator,
cached_unary_operator)
from ..constraints import (AffineConstraint, ComplexAffineConstraint,
ComplexLMIConstraint, Constraint, LMIConstraint)
from ..formatting import detect_range
from ..legacy import deprecated
from .data import (blend_shapes, convert_operands, cvx2np, cvxopt_vcat,
left_kronecker_I, load_data, load_dense_data, load_shape,
right_kronecker_I, sparse_quadruple)
from .expression import (Expression, ExpressionType, refine_operands,
validate_prediction)
_API_START = api_start(globals())
# -------------------------------
RELATIVE_HERMITIANNESS_TOLERANCE = 1e-10
r"""Relative tolerance when checking if an expression is hermitian.
A square (complex) affine expression :math:`A` is considered hermitian if
.. math::
\max_{1 \leq i, j \leq n} |(A - A^H)_{ij}|
\leq
\varepsilon \max_{1 \leq i, j \leq n} |A_{ij}|
where :math:`\varepsilon` is this tolerance.
Used by :meth:`ComplexAffineExpression.hermitian`.
"""
[docs]class ComplexAffineExpression(Expression):
"""A multidimensional complex affine expression.
Base class for the real :class:`AffineExpression`.
"""
# --------------------------------------------------------------------------
# Initialization and factory methods.
# --------------------------------------------------------------------------
[docs] def __init__(self, string, shape=(1, 1), coefficients={}, constant=0j):
"""Initialize a (complex) affine expression.
This constructor is meant for internal use. As a user, you will most
likely want to build expressions starting from
:mod:`~picos.expressions.variables` or a :func:`~picos.Constant`.
:param str string: A symbolic string description.
:param shape: Shape of a vector or matrix expression.
:type shape: int or tuple or list
:param dict coefficients: The linear part of the affine expression; maps
the vectorization of :class:`variables <.variables.BaseVariable>` to
their coefficients.
:param constant: The constant part of the affine expression.
.. warning::
If the given coefficients and constant are already of the desired
(numeric) type and shape, they are stored by reference. Modifying
such data can lead to unexpected results as PICOS expressions are
supposed to be immutable (to allow
`caching <https://en.wikipedia.org/wiki/Cache_(computing)>`_ of
results).
If you create affine expressions by any other means than this
constructor, PICOS makes a copy of your data to prevent future
modifications to it from causing inconsistencies.
"""
from .variables import BaseVariable
shape = load_shape(shape)
length = shape[0]*shape[1]
if not isinstance(coefficients, dict):
raise TypeError("Coefficients of {} must be given as a dict."
.format(type(self).__name__))
# Store shape.
self._shape = shape
# Store coefficients.
self._coefs = {}
for var, coef in coefficients.items():
# Do not store coefficients of zero.
if not coef:
continue
if not isinstance(var, BaseVariable):
raise TypeError("Coefficients of {} must be indexed by "
"variables, not objects of type {}.".format(
type(self).__name__, type(var).__name__))
self._coefs[var] = load_data(
coef, (length, var.dim), self._typecode, alwaysCopy=False)[0]
# Store vectorized constant.
self._const = load_data(
constant, length, self._typecode, alwaysCopy=False)[0]
# Determine the type string.
typeBaseStr = self._get_type_string_base()
if "{}" in typeBaseStr:
if self._coefs:
if self._const:
typeBaseStr = typeBaseStr.format("Affine Expression")
else:
typeBaseStr = typeBaseStr.format("Linear Expression")
else:
typeBaseStr = typeBaseStr.format("Constant")
typeStr = "{} {}".format(glyphs.shape(shape), typeBaseStr)
Expression.__init__(self, typeStr, string)
@classmethod
def _get_basetype(cls):
"""Return the first order :class:`Expression` subclass this is based on.
Enables subclass objects (such as variables) to behave like the returned
type with respect to algebraic operations. For instance, the sum of
two :class:`ComplexVariable` is a :class:`ComplexAffineExpression`.
"""
return ComplexAffineExpression
@classmethod
def _get_typecode(cls):
"""See :meth:`_get_basetype`."""
return "z"
@classmethod
def _get_type_string_base(cls):
"""See :meth:`_get_basetype`."""
return "Complex {}"
_basetype = property(lambda self: self._get_basetype())
_typecode = property(lambda self: self._get_typecode())
_typeStrBase = property(lambda self: self._get_type_string_base())
[docs] @classmethod
def from_constant(cls, constant, shape=None, name=None):
"""Create a class instance from the given numeric constant.
Loads the given constant as a PICOS expression, optionally broadcasted
or reshaped to the given shape and named as specified.
See :func:`~.data.load_data` for supported data formats and broadcasting
and reshaping rules.
Unlike :func:`Constant`, this class method always creates an instance of
the class that it is called on, instead of tailoring towards the numeric
type of the data.
.. note::
When an operation involves both a PICOS expression and a constant
value of another type, PICOS converts the constant on the fly so
that you rarely need to use this method.
"""
constant, string = load_data(constant, shape, cls._get_typecode())
return cls(name if name else string, constant.size, {}, constant)
[docs] @classmethod
def zero(cls, shape=(1, 1)):
"""Return a constant zero expression of given shape."""
shape = load_shape(shape)
string = glyphs.scalar(0) if shape == (1, 1) else glyphs.matrix(0)
return cls(string, shape)
# --------------------------------------------------------------------------
# Abstract method implementations and method overridings, except _predict.
# --------------------------------------------------------------------------
@cached_unary_operator
def _get_refined(self):
if self.isreal:
return AffineExpression(
self._symbStr, self._shape, self._coefs, self._const)
else:
return self
Subtype = namedtuple("Subtype", ("shape", "constant", "nonneg"))
Subtype.dim = property(lambda self: self.shape[0] * self.shape[1])
def _get_subtype(self):
nonneg = self.constant and self.isreal \
and all(x >= 0 for x in self.value_as_matrix)
return self.Subtype(self._shape, self.constant, nonneg)
def _get_value(self):
# Create a vectorized copy of the constant term.
value = self._const[:]
for var, coef in self._coefs.items():
summand = coef * var._get_internal_value()
if type(value) == type(summand):
value += summand
else:
# Exactly one of the matrices is sparse.
value = value + summand
# Resize the value to the proper shape.
value.size = self.shape
return value
def _devectorize(self, value):
value.size = self._shape
return value
@convert_operands(sameShape=True, allowNone=True)
def _set_value(self, value):
if value is None:
for var in self._coefs:
var.value = None
return
# Since all variables are real-valued, prevent NumPy from finding
# complex-valued solutions that do not actually work.
(self.real // self.imag).renamed(self.string).value \
= (value.real // value.imag)
def _get_shape(self):
return self._shape
@cached_unary_operator
def _get_variables(self):
return frozenset(var for var in self._coefs)
def _is_convex(self):
return True
def _is_concave(self):
return True
def _replace_variables(self, var_map):
# Handle the base case where the affine expression is a variable.
if self in var_map:
return var_map[self]
name_map = {old.name: new.name for old, new in var_map.items()}
string = self.string
if isinstance(string, glyphs.GlStr):
string = string.reglyphed(name_map)
elif string in name_map:
# Handle corner cases like x + 0 for a variable x.
string = name_map[string]
coefs = {var_map[var]: coef for var, coef in self._coefs.items()}
return self._basetype(string, self._shape, coefs, self._const)
# --------------------------------------------------------------------------
# Python special method implementations, except constraint-creating ones.
# --------------------------------------------------------------------------
def __len__(self):
# Faster version that overrides Expression.__len__.
return self._shape[0] * self._shape[1]
def __getitem__(self, index):
def slice2range(s, length):
"""Transform a :class:`slice` to a :class:`range`."""
assert isinstance(s, slice)
# Plug in slice's default values.
ss = s.step if s.step else 1
if ss > 0:
sa = s.start if s.start is not None else 0
sb = s.stop if s.stop is not None else length
else:
assert ss < 0
sa = s.start if s.start is not None else length - 1
sb = s.stop # Keep it None as -1 would mean length - 1.
# Wrap negative indices (once).
ra = length + sa if sa < 0 else sa
if sb is None:
# This is the only case where we give a negative index to range.
rb = -1
else:
rb = length + sb if sb < 0 else sb
# Clamp out-of-bound indices.
ra = min(max(0, ra), length - 1)
rb = min(max(-1, rb), length)
r = range(ra, rb, ss)
if not r:
raise IndexError("Empty slice.")
return r
def range2slice(r, length):
"""Transform a :class:`range` to a :class:`slice`, if possible.
:raises ValueError: If the input cannot be expressed as a slice.
"""
assert isinstance(r, range)
if not r:
raise IndexError("Empty range.")
ra = r.start
rb = r.stop
rs = r.step
if rs > 0:
if ra < 0 or rb > length:
raise ValueError(
"Out-of-bounds range cannot be represented as a slice.")
else:
assert rs < 0
if ra >= length or rb < -1:
raise ValueError(
"Out-of-bounds range cannot be represented as a slice.")
if rb == -1:
rb = None
return slice(ra, rb, rs)
def list2slice(l, length):
"""Transform a :class:`list` to a :class:`slice`, if possible.
:raises TypeError: If the input is not an integer sequence.
:raises ValueError: If the input cannot be expressed as a slice.
"""
r = detect_range(l)
if isinstance(r, list): # Python 2
raise RuntimeError(
"Unable to convert from list to slice because the "
"intermediate range object does not exist under Python 2.")
return range2slice(r, length)
def slice2str(s):
"""Return the short string that produced a :class:`slice`."""
assert isinstance(s, slice)
startStr = str(s.start) if s.start is not None else ""
stopStr = str(s.stop) if s.stop is not None else ""
if s.step in (None, 1):
if s.start is not None and s.stop is not None \
and s.stop == s.start + 1:
return startStr
else:
return "{}:{}".format(startStr, stopStr)
else:
return "{}:{}:{}".format(startStr, stopStr, str(s.step))
def list2str(l, length):
"""Return a short string represnetation of a :class:`list`."""
assert isinstance(l, list)
# Extract integers wrapped in a list.
if len(l) == 1:
return str(l[0])
# Try to convert the list to a slice.
try:
l = list2slice(l, length)
except (ValueError, RuntimeError):
pass
if isinstance(l, list):
if len(l) > 4:
return glyphs.shortint(l[0], l[-1])
else:
return str(l).replace(" ", "")
else:
return slice2str(l)
def any2str(a, length):
if isinstance(a, slice):
return slice2str(a)
elif isinstance(a, list):
return list2str(a, length)
else:
assert False
m, n = self._shape
indexStr = None
isIntList = False
# Turn the index expression into a mutable list of one index per axis.
if isinstance(index, tuple): # Multiple axis slicing.
index = list(index)
elif isinstance(index, dict): # Arbitrary matrix element selection.
if len(index) != 2:
"If slicing with a dictionary, there must be exactly two keys."
keys = index.keys()
i, j = min(keys), max(keys)
I, J = index[i], index[j]
try:
I = load_dense_data(I, typecode="i", alwaysCopy=False)[0]
J = load_dense_data(J, typecode="i", alwaysCopy=False)[0]
if 1 not in I.size or 1 not in J.size:
raise TypeError("At least one of the objects is not flat "
"but a proper matrix.")
if len(I) != len(J):
raise TypeError("The objects do not have the same length.")
I, J = list(I), list(J)
except Exception as error:
raise TypeError("Loading a sparse index vector pair for {} from"
" objects of type {} and {} failed: {}".format(
self.string, type(I).__name__, type(J).__name__, error))
# Represent the selection as a global index list.
index = [[i + j*m for i, j in zip(I, J)]]
# Use a special index string.
indexStr = glyphs.size(list2str(I, n), list2str(J, m))
# Don't invoke load_dense_data on the list.
isIntList = True
else: # Global indexing.
index = [index]
# Make sure either global or row/column indexing is used.
if not index:
raise IndexError("Empty index.")
elif len(index) > 2:
raise IndexError(
"PICOS expressions do not have a third axis to slice.")
# Turn the indices for each axis into either a slice or a list.
for axis, selection in enumerate(index):
# Convert anything that is not a slice, including scalars and lists
# that are not confirmed integer, to an integer row or column
# vector, then (back) to a list.
if not isIntList and not isinstance(selection, slice):
try:
matrix = load_dense_data(
selection, typecode="i", alwaysCopy=False)[0]
if 1 not in matrix.size:
raise TypeError("The object is not flat but a {} shaped"
" matrix.".format(glyphs.shape(matrix.size)))
selection = list(matrix)
except Exception as error:
raise TypeError("Loading a slicing index vector for axis {}"
" of {} from an object of type {} failed: {}".format(
axis, self.string, type(selection).__name__, error))
index[axis] = selection
# Build index string, retrieve new shape, finalize index.
if len(index) == 1: # Global indexing.
index = index[0]
if isinstance(index, slice):
shape = len(slice2range(index, len(self)))
else:
shape = len(index)
if indexStr is None:
indexStr = any2str(index, len(self))
else: # Multiple axis slicing.
if indexStr is None:
indexStr = "{},{}".format(
any2str(index[0], m), any2str(index[1], n))
# Convert to a global index list.
RC, shape = [], []
for axis, selection in enumerate(index):
k = self._shape[axis]
if isinstance(selection, slice):
# Turn the slice into an iterable range.
selection = slice2range(selection, k)
# All indices from a slice are nonnegative.
assert all(i >= 0 for i in selection)
if isinstance(selection, list):
# Wrap once. This is consistent with CVXOPT.
selection = [i if i >= 0 else k + i for i in selection]
# Perform a partial out-of-bounds check.
if any(i < 0 for i in selection):
raise IndexError(
"Out-of-bounds access along axis {}.".format(axis))
# Complete the check for out-of-bounds access.
if any(i >= k for i in selection):
raise IndexError(
"Out-of-bounds access along axis {}.".format(axis))
RC.append(selection)
shape.append(len(selection))
rows, cols = RC
index = [i + j*m for j in cols for i in rows]
# Finalize the string.
string = glyphs.slice(self.string, indexStr)
# Retrieve new coefficients and constant term.
try:
coefs = {var: coef[index, :] for var, coef in self._coefs.items()}
const = self._const[index]
except IndexError as error:
raise IndexError("Invalid expression slice {}: {}".format(
string, error))
return self._basetype(string, shape, coefs, const)
def _common_basetype(self, other):
assert isinstance(other, ComplexAffineExpression)
if issubclass(other._basetype, self._basetype):
return self._basetype
else:
return other._basetype
@convert_operands(sameShape=True)
@refine_operands(stop_at_affine=True)
def __add__(self, other):
"""Denote addition from the right hand side."""
if not isinstance(other, ComplexAffineExpression):
return NotImplemented
string = glyphs.clever_add(self.string, other.string)
coefs = {}
for var, coef in self._coefs.items():
coefs[var] = coef+other._coefs[var] if var in other._coefs else coef
for var, coef in other._coefs.items():
coefs.setdefault(var, coef)
const = self._const + other._const
return self._common_basetype(other)(string, self._shape, coefs, const)
@convert_operands(sameShape=True)
@refine_operands(stop_at_affine=True)
def __radd__(self, other):
"""Denote addition from the left hand side."""
if not isinstance(other, ComplexAffineExpression):
return NotImplemented
return other.__add__(self)
@convert_operands(sameShape=True)
@refine_operands(stop_at_affine=True)
def __sub__(self, other):
"""Denote substraction from the right hand side."""
if not isinstance(other, ComplexAffineExpression):
return NotImplemented
string = glyphs.clever_sub(self.string, other.string)
coefs = {}
for var, coef in self._coefs.items():
coefs[var] = coef-other._coefs[var] if var in other._coefs else coef
for var, coef in other._coefs.items():
coefs.setdefault(var, -coef)
const = self._const - other._const
return self._common_basetype(other)(string, self._shape, coefs, const)
@convert_operands(sameShape=True)
@refine_operands(stop_at_affine=True)
def __rsub__(self, other):
"""Denote substraction with self on the right hand side."""
if not isinstance(other, ComplexAffineExpression):
return NotImplemented
return other.__sub__(self)
@cached_selfinverse_unary_operator
def __neg__(self):
"""Denote negation."""
string = glyphs.clever_neg(self.string)
coefs = {var: -coef for var, coef in self._coefs.items()}
const = -self._const
return self._basetype(string, self._shape, coefs, const)
@convert_operands(sameShape=True)
@refine_operands(stop_at_affine=True)
def __or__(self, other):
r"""Denote the scalar product with self on the left hand side.
For (complex) vectors :math:`a` and :math:`b` this is the dot product
.. math::
(a \mid b)
&= \langle a, b \rangle \\
&= a \cdot b \\
&= b^H a.
For (complex) matrices :math:`A` and :math:`B` this is the Frobenius
inner product
.. math::
(A \mid B)
&= \langle A, B \rangle_F \\
&= A : B \\
&= \operatorname{tr}(B^H A) \\
&= \operatorname{vec}(B)^H \operatorname{vec}(\overline{A})
.. note::
Write ``(A|B)`` instead of ``A|B`` for the scalar product of ``A``
and ``B`` to obtain correct operator binding within a larger
expression context.
"""
from .exp_quadratic import QuadraticExpression
if not isinstance(other, ComplexAffineExpression):
return NotImplemented
# Make a string description.
string = glyphs.clever_dotp(
self.string, other.string, other.complex, self.scalar)
# Handle the complex case: Conjugate the right hand side.
other = other.conj
if other.constant:
factor = other._const.T
coefs = {var: factor * coef for var, coef in self._coefs.items()}
const = factor * self._const
elif self.constant:
factor = self._const.T
coefs = {var: factor * coef for var, coef in other._coefs.items()}
const = factor * other._const
else:
# Compute the affine part of the product.
affString = glyphs.affpart(string)
affCoefs = {}
for var in set(self._coefs.keys()).union(other._coefs.keys()):
if var not in other._coefs:
affCoefs[var] = other._const.T * self._coefs[var]
elif var not in self._coefs:
affCoefs[var] = self._const.T * other._coefs[var]
else:
affCoefs[var] = other._const.T * self._coefs[var] + \
self._const.T * other._coefs[var]
affConst = self._const.T * other._const
affPart = self._common_basetype(other)(
affString, (1, 1), affCoefs, affConst)
# Compute the quadratic part of the product.
quadPart = {(v, w): self._coefs[v].T * other._coefs[w]
for v in self._coefs for w in other._coefs}
# Don't create quadratic expressions without a quadratic part.
if not any(quadPart.values()):
affPart._symbStr = string
return affPart
# Remember a factorization into two real scalars if applicable.
# NOTE: If the user enters a multiplication a*b of two scalar affine
# expressions, then we have, at this point, self == a.T == a
# and other == b.conj.conj == b.
if len(self) == 1 and len(other) == 1 \
and self.isreal and other.isreal:
factors = (self.refined, other.refined)
else:
factors = None
return QuadraticExpression(
string, quadPart, affPart, scalarFactors=factors)
return self._common_basetype(other)(string, (1, 1), coefs, const)
@convert_operands(sameShape=True)
@refine_operands(stop_at_affine=True)
def __ror__(self, other):
"""Denote the scalar product with self on the right hand side."""
if not isinstance(other, ComplexAffineExpression):
return NotImplemented
return other.__or__(self)
@convert_operands(rMatMul=True)
@refine_operands(stop_at_affine=True)
def __mul__(self, other):
"""Denote matrix multiplication from the right hand side."""
if not isinstance(other, ComplexAffineExpression):
return NotImplemented
string = glyphs.clever_mul(self.string, other.string)
shape = None
if other.constant:
if other._shape == (1, 1):
factor = other._const[0]
shape = self._shape
elif self._shape == (1, 1):
factor = other._const
shape = other._shape
else:
# Use the identity vec(AB) = (Bᵀ ⊗ I)vec(A).
factor = right_kronecker_I(other._const, self._shape[0],
reshape=other._shape, postT=True)
coefs = {var: factor * coef for var, coef in self._coefs.items()}
const = factor * self._const
elif self.constant:
if self._shape == (1, 1):
factor = self._const[0]
shape = other._shape
elif other._shape == (1, 1):
factor = self._const
shape = self._shape
else:
# Use the identity vec(AB) = (I ⊗ A)vec(B).
factor = left_kronecker_I(
self._const, other._shape[1], reshape=self._shape)
coefs = {var: factor * coef for var, coef in other._coefs.items()}
const = factor * other._const
else:
# If the result is scalar, allow for quadratic terms.
if self._shape[0] == 1 and other._shape[1] == 1 \
and self._shape[1] == other._shape[0]:
result = self.T.__or__(other.conj)
# NOTE: __or__ always creates a fresh expression.
result._symbStr = string
return result
else:
raise NotImplementedError(
"PICOS does not support multidimensional quadratic "
"expressions at this point. More precisely, one factor must"
" be constant or the result must be scalar.")
if shape is None:
shape = (self._shape[0], other._shape[1])
return self._common_basetype(other)(string, shape, coefs, const)
@convert_operands(lMatMul=True)
@refine_operands(stop_at_affine=True)
def __rmul__(self, other):
"""Denote matrix multiplication from the left hand side."""
if not isinstance(other, ComplexAffineExpression):
return NotImplemented
return other.__mul__(self)
@convert_operands(sameShape=True)
@refine_operands(stop_at_affine=True)
def __xor__(self, other):
"""Denote the elementwise (or Hadamard) product."""
if not isinstance(other, ComplexAffineExpression):
return NotImplemented
string = glyphs.hadamard(self.string, other.string)
if other.constant:
factor = cvxopt.spdiag(other._const)
coefs = {var: factor * coef for var, coef in self._coefs.items()}
const = factor * self._const
elif self.constant:
factor = cvxopt.spdiag(self._const)
coefs = {var: factor * coef for var, coef in other._coefs.items()}
const = factor * other._const
else:
# If the result is scalar, allow for quadratic terms.
if self._shape == (1, 1):
result = self.__or__(other.conj)
# NOTE: __or__ always creates a fresh expression.
result._symbStr = string
return result
else:
raise NotImplementedError(
"PICOS does not support multidimensional quadratic "
"expressions at this point. More precisely, one factor must"
" be constant or the result must be scalar.")
return self._common_basetype(other)(string, self._shape, coefs, const)
@convert_operands(sameShape=True)
@refine_operands(stop_at_affine=True)
def __rxor__(self, other):
"""Denote the elementwise (or Hadamard) product."""
if not isinstance(other, ComplexAffineExpression):
return NotImplemented
return other.__xor__(self)
@convert_operands()
@refine_operands(stop_at_affine=True)
def __matmul__(self, other):
"""Denote the Kronecker product from the right hand side."""
def vectorized_kron(A, B):
m, n = self._shape
p, q = other._shape
mn, r = A.size
pq, s = B.size
mnpq = mn*pq
assert 1 in (r, s)
assert m*n == mn
assert p*q == pq
if r < s:
A, P = cvx2np(A, (m, n)), []
for j in range(s):
Bj = cvx2np(B[:, j], (p, q))
P.append(numpy.reshape(numpy.kron(A, Bj), (mnpq, 1), "F"))
else:
B, P = cvx2np(B, (p, q)), []
for j in range(r):
Aj = cvx2np(A[:, j], (m, n))
P.append(numpy.reshape(numpy.kron(Aj, B), (mnpq, 1), "F"))
assert P
return load_data(numpy.vstack(P), typecode=tc)[0]
if not isinstance(other, ComplexAffineExpression):
return NotImplemented
cls = self._common_basetype(other)
tc = cls._get_typecode()
string = glyphs.kron(self.string, other.string)
m, n = self._shape
p, q = other._shape
shape = (m*p, n*q)
if other.constant:
coefs = {var: vectorized_kron(coef, other._const)
for var, coef in self._coefs.items()}
elif self.constant:
coefs = {var: vectorized_kron(self._const, coef)
for var, coef in other._coefs.items()}
else:
raise NotImplementedError(
"PICOS does not support multidimensional quadratic expressions "
"at this point. More precisely, at least one factor of the "
"Kronecker product must be constant.")
const = vectorized_kron(self._const, other._const)
return cls(string, shape, coefs, const)
@convert_operands()
@refine_operands(stop_at_affine=True)
def __rmatmul__(self, other):
"""Denote the Kronecker product from the left hand side."""
if not isinstance(other, ComplexAffineExpression):
return NotImplemented
return other.__matmul__(self)
[docs] def kron(self, other):
"""Denote the Kronecker product from the right hand side.
Python 3 users can use the infix ``@`` operator instead.
"""
# Python 2 does not have the @ operator.
return self.__matmul__(other)
[docs] def leftkron(self, other):
"""Denote the Kronecker product from the left hand side.
Python 3 users can use the infix ``@`` operator instead.
"""
# Python 2 does not have the @ operator.
return self.__rmatmul__(other)
@convert_operands(scalarRHS=True)
@refine_operands(stop_at_affine=True)
def __truediv__(self, other):
"""Denote elementwise division by a scalar constant."""
if not isinstance(other, ComplexAffineExpression):
return NotImplemented
if not other.constant:
raise TypeError(
"You may only divide {} by a constant.".format(self.string))
if other.is0:
raise ZeroDivisionError(
"Tried to divide {} by zero.".format(self.string))
divisor = other._const[0]
string = glyphs.div(self.string, other.string)
coefs = {var: coef / divisor for var, coef in self._coefs.items()}
const = self._const / divisor
return self._common_basetype(other)(string, self._shape, coefs, const)
@convert_operands(scalarLHS=True)
@refine_operands(stop_at_affine=True)
def __rtruediv__(self, other):
"""Denote elementwise division by scalar constant self."""
if not isinstance(other, ComplexAffineExpression):
return NotImplemented
return other.__div__(self)
def __div__(self, other):
"""Denote elementwise division by a scalar constant."""
# Python 2 uses this instead of __truediv__.
return self.__truediv__(other)
def __rdiv__(self, other):
"""Denote elementwise division by scalar constant self."""
# Python 2 uses this instead of __rtruediv__.
return self.__rtruediv__(other)
@convert_operands(horiCat=True)
@refine_operands(stop_at_affine=True)
def __and__(self, other):
"""Denote horizontal concatenation from the right hand side."""
if not isinstance(other, ComplexAffineExpression):
return NotImplemented
string = glyphs.matrix_cat(self.string, other.string, horizontal=True)
shape = (self._shape[0], self._shape[1] + other._shape[1])
coefs = {}
for var in set(self._coefs.keys()).union(other._coefs.keys()):
lhs = self._coefs[var] if var in self._coefs else cvxopt.spmatrix(
[], [], [], (len(self), other._coefs[var].size[1]))
rhs = other._coefs[var] if var in other._coefs else cvxopt.spmatrix(
[], [], [], (len(other), self._coefs[var].size[1]))
coefs[var] = cvxopt_vcat([lhs, rhs])
const = cvxopt_vcat([self._const, other._const])
return self._common_basetype(other)(string, shape, coefs, const)
@convert_operands(horiCat=True)
@refine_operands(stop_at_affine=True)
def __rand__(self, other):
"""Denote horizontal concatenation from the left hand side."""
if not isinstance(other, ComplexAffineExpression):
return NotImplemented
return other.__and__(self)
@convert_operands(vertCat=True)
@refine_operands(stop_at_affine=True)
def __floordiv__(self, other):
"""Denote vertical concatenation from below."""
def interleave_columns(upper, lower, upperRows, lowerRows, cols):
p, q = upperRows, lowerRows
return [column for columnPairs in [
(upper[j*p:j*p+p, :], lower[j*q:j*q+q, :]) for j in range(cols)]
for column in columnPairs]
if not isinstance(other, ComplexAffineExpression):
return NotImplemented
string = glyphs.matrix_cat(self.string, other.string, horizontal=False)
p, q, c = self._shape[0], other._shape[0], self._shape[1]
shape = (p + q, c)
coefs = {}
for var in set(self._coefs.keys()).union(other._coefs.keys()):
upr = self._coefs[var] if var in self._coefs else cvxopt.spmatrix(
[], [], [], (len(self), other._coefs[var].size[1]))
lwr = other._coefs[var] if var in other._coefs else cvxopt.spmatrix(
[], [], [], (len(other), self._coefs[var].size[1]))
coefs[var] = cvxopt_vcat(interleave_columns(upr, lwr, p, q, c))
const = cvxopt_vcat(
interleave_columns(self._const, other._const, p, q, c))
return self._common_basetype(other)(string, shape, coefs, const)
@convert_operands(vertCat=True)
@refine_operands(stop_at_affine=True)
def __rfloordiv__(self, other):
"""Denote vertical concatenation from above."""
if not isinstance(other, ComplexAffineExpression):
return NotImplemented
return other.__floordiv__(self)
@convert_operands(scalarRHS=True)
@refine_operands() # Refine both sides to real if possible.
def __pow__(self, other):
"""Denote exponentiation."""
from .exp_powtrace import PowerTrace
if not isinstance(other, AffineExpression):
return NotImplemented
if not self.scalar:
raise TypeError("May only exponentiate a scalar expression.")
if not other.constant:
raise TypeError("The exponent must be constant.")
exponent = other.value
if exponent == 2:
return (self | self) # Works for complex base.
else:
return PowerTrace(self, exponent) # Errors on complex base.
@cached_unary_operator
def __abs__(self):
"""Denote the Euclidean or Frobenius norm of the expression."""
from . import Norm
return Norm(self)
# --------------------------------------------------------------------------
# Properties and functions that describe the expression.
# --------------------------------------------------------------------------
@cached_property
def _sparse_coefs(self):
"""Linear coefficients (re-)cast to sparse matrices."""
return {var: cvxopt.sparse(coef) for var, coef in self._coefs.items()}
@property
def scalar(self):
"""Whether this is a scalar affine expression."""
return self._shape == (1, 1)
@property
def constant(self):
"""Whether the expression involves no variables."""
return not self._coefs
@property
def square(self):
"""Whether the expression is a square matrix."""
return self._shape[0] == self._shape[1]
@cached_property
def hermitian(self): # noqa (D402 thinks this includes a signature)
"""Whether the expression is a hermitian (or symmetric) matrix.
A tolerance is used to compensate for small numeric errors that can
occur when creating hermtian expressions by means of matrix
multiplication, see
:data:`~.exp_affine.RELATIVE_HERMITIANNESS_TOLERANCE`.
If PICOS still rejects your expression as not hermitian (or as not
symmetric), you can use :meth:`hermitianized` to correct larger numeric
errors or the effects of noisy data.
"""
return self.equals(self.H, relTol=RELATIVE_HERMITIANNESS_TOLERANCE)
@cached_property
def is0(self):
"""Whether this is a constant scalar, vector or matrix of all zeros."""
# Must be constant.
if self._coefs:
return False
# May not have a nonzero constant term.
if self._const:
return False
return True
@cached_property
def is1(self):
"""Whether this is a constant scalar or vector of all ones."""
# Must be a scalar or vector.
if self._shape[0] != 1 and self._shape[1] != 1:
return False
# Must be constant.
if self._coefs:
return False
# Constant term must be all ones.
for index in range(len(self)):
if self._const[index] != 1:
return False
return True
@cached_property
def isI(self):
"""Whether this is a constant identity matrix."""
m, n = self._shape
# Must be a square matrix.
if m != n:
return False
# Must be constant.
if self._coefs:
return False
# Constant term must be the identity.
for row in range(m):
for col in range(n):
if self._const[col*m + row] != int(row == col):
return False
return True
@cached_property
def isreal(self):
"""Whether the expression is always real-valued."""
if any(self._const.imag()):
return False
for coef in self._coefs.values():
if any(coef.imag()):
return False
return True
@property
def complex(self):
"""Whether the expression can be complex-valued."""
return not self.isreal
[docs] @convert_operands()
def equals(self, other, absTol=None, relTol=None):
"""Check mathematical equality with another affine expression.
Checks whether the affine expression involves the same variables with
the same coefficients and an equal constant term as the given PICOS
(affine) expression or constant term.
The type of both expressions may differ. In particular, a
:class:`ComplexAffineExpression` with real coefficients and constant can
be equal to an :class:`AffineExpression`.
If the argument is a constant term, no reshaping or broadcasting is used
to bring it to the same shape as this expression. In particular,
- ``0`` refers to a scalar zero (see also :meth:`is0`),
- lists and tuples are treated as column vectors and
- algebraic strings must specify a shape (see
:func:`~.data.load_data`).
:param other: Another PICOS expression or a constant numeric data value
supported by :func:`~.data.load_data`.
:param absTol: As long as all absolute differences between scalar
entries of the coefficient matrices and the constant terms being
compared does not exceed this bound, consider the expressions equal.
:param relTol: As long as all absolute differences between scalar
entries of the coefficient matrices and the constant terms being
compared divided by the maximum absolute value found in either term
does not exceed this bound, consider the expressions equal.
:Example:
>>> from picos import Constant
>>> A = Constant("A", 0, (5,5))
>>> repr(A)
'<5×5 Real Constant: A>'
>>> A.is0
True
>>> A.equals(0)
False
>>> A.equals("|0|(5,5)")
True
>>> repr(A*1j)
'<5×5 Complex Constant: A·1j>'
>>> A.equals(A*1j)
True
"""
def differ(A, B):
Z = A - B
if Z:
if not absTol and not relTol:
return True
M = max(abs(Z))
if relTol:
N = max(max(abs(A)), max(abs(B)))
if absTol and relTol:
if M > absTol and M / N > relTol:
return True
elif absTol:
if M > absTol:
return True
else:
if M / N > relTol:
return True
return False
if self is other:
return True
if not isinstance(other, ComplexAffineExpression):
return False
if self._shape != other._shape:
return False
if differ(self._const, other._const):
return False
for var in other._coefs:
if var not in self._coefs:
return False
for var, coef in self._coefs.items():
if var not in other._coefs:
return False
if differ(coef, other._coefs[var]):
return False
return True
# --------------------------------------------------------------------------
# Methods and properties that return modified copies.
# --------------------------------------------------------------------------
[docs] def renamed(self, string):
"""Return the expression with a modified string description."""
return self._basetype(string, self._shape, self._coefs, self._const)
[docs] def reshaped(self, shape):
"""Return the expression reshaped in column-major order.
:Example:
>>> from picos import Constant
>>> C = Constant("C", range(6), (2, 3))
>>> print(C)
[ 0.00e+00 2.00e+00 4.00e+00]
[ 1.00e+00 3.00e+00 5.00e+00]
>>> print(C.reshaped((3, 2)))
[ 0.00e+00 3.00e+00]
[ 1.00e+00 4.00e+00]
[ 2.00e+00 5.00e+00]
"""
shape = load_shape(shape, wildcards=True)
if shape == self._shape:
return self
elif shape == (None, None):
return self
length = len(self)
if shape[0] is None:
shape = (length // shape[1], shape[1])
elif shape[1] is None:
shape = (shape[0], length // shape[0])
if shape[0]*shape[1] != length:
raise ValueError("Can only reshape to a matrix of same size.")
string = glyphs.reshaped(self.string, glyphs.shape(shape))
return self._basetype(string, shape, self._coefs, self._const)
[docs] def broadcasted(self, shape):
"""Return the expression broadcasted to the given shape.
:Example:
>>> from picos import Constant
>>> C = Constant("C", range(6), (2, 3))
>>> print(C)
[ 0.00e+00 2.00e+00 4.00e+00]
[ 1.00e+00 3.00e+00 5.00e+00]
>>> print(C.broadcasted((6, 6)))
[ 0.00e+00 2.00e+00 4.00e+00 0.00e+00 2.00e+00 4.00e+00]
[ 1.00e+00 3.00e+00 5.00e+00 1.00e+00 3.00e+00 5.00e+00]
[ 0.00e+00 2.00e+00 4.00e+00 0.00e+00 2.00e+00 4.00e+00]
[ 1.00e+00 3.00e+00 5.00e+00 1.00e+00 3.00e+00 5.00e+00]
[ 0.00e+00 2.00e+00 4.00e+00 0.00e+00 2.00e+00 4.00e+00]
[ 1.00e+00 3.00e+00 5.00e+00 1.00e+00 3.00e+00 5.00e+00]
"""
shape = load_shape(shape, wildcards=True)
shape = blend_shapes(shape, self._shape)
if shape == self._shape:
return self
vdup = shape[0] // self._shape[0]
hdup = shape[1] // self._shape[1]
if (self._shape[0] * vdup, self._shape[1] * hdup) != shape:
raise ValueError("Cannot broadcast from shape {} to {}."
.format(glyphs.shape(self._shape), glyphs.shape(shape)))
if self._shape == (1, 1):
string = glyphs.matrix(self.string)
return (self * cvxopt.matrix(1.0, shape)).renamed(string)
string = glyphs.bcasted(self.string, glyphs.shape(shape))
return self.leftkron(cvxopt.matrix(1.0, (vdup, hdup))).renamed(string)
[docs] def reshaped_or_broadcasted(self, shape):
"""Return the expression :meth:`reshaped` or :meth:`broadcasted`.
Unlike with :meth:`reshaped` and :meth:`broadcasted`, the target shape
may not contain a wildcard character.
If the wildcard-free target shape has the same number of elements as
the current shape, then this is the same as :meth:`reshaped`, otherwise
it is the same as :meth:`broadcasted`.
"""
shape = load_shape(shape, wildcards=False)
try:
if shape[0]*shape[1] == len(self):
return self.reshaped(shape)
else:
return self.broadcasted(shape)
except ValueError:
raise ValueError("Cannot reshape or broadcast from shape {} to {}."
.format(glyphs.shape(self._shape), glyphs.shape(shape)))
@cached_property
def hermitianized(self):
r"""The expression projected onto the subspace of hermitian matrices.
For a square (complex) affine expression :math:`A`, this is
:math:`\frac{1}{2}(A + A^H)`.
If the expression is not complex, then this is the projection onto the
subspace of symmetric matrices.
"""
if not self.square:
raise TypeError("Cannot hermitianize non-square {}.".format(self))
return (self + self.H)/2
@cached_property
def real(self):
"""Real part of a complex affine expression."""
return AffineExpression(
glyphs.real(self.string),
self._shape,
{var: coef.real() for var, coef in self._coefs.items()},
self._const.real())
@cached_property
def imag(self):
"""Imaginary part of a complex affine expression."""
return AffineExpression(
glyphs.imag(self.string),
self._shape,
{var: coef.imag() for var, coef in self._coefs.items()},
self._const.imag())
@cached_property
def lin(self):
"""Linear part of the affine expression."""
return self._basetype(
glyphs.linpart(self._symbStr), self._shape, self._coefs)
@cached_property
def cst(self):
"""Constant part of the affine expression."""
return self._basetype(
glyphs.cstpart(self._symbStr), self._shape, constant=self._const)
@cached_selfinverse_property
def T(self):
"""Matrix transpose."""
if len(self) == 1:
return self
# Construct a commutation matrix to transpose the vectorized data.
n, m = self._shape
d = n * m
V = [1]*d
I = range(d)
J = [(k % m)*n + k // m for k in I]
K = cvxopt.spmatrix(V, I, J, (d, d), self._typecode)
string = glyphs.transp(self.string)
shape = (self._shape[1], self._shape[0])
coefs = {var: K * coef for var, coef in self._coefs.items()}
const = K * self._const
return self._basetype(string, shape, coefs, const)
@cached_selfinverse_property
def conj(self):
"""Complex conjugate."""
string = glyphs.conj(self.string)
coefs = {var: coef.H.T for var, coef in self._coefs.items()}
const = self._const.H.T
return self._basetype(string, self._shape, coefs, const)
@cached_selfinverse_property
def H(self):
"""Conjugate (or Hermitian) transpose."""
return self.T.conj.renamed(glyphs.htransp(self._symbStr))
def _square_equal_subsystem_dims(self, diagLen):
"""Support :func:`partial_trace` and :func:`partial_transpose`."""
m, n = self._shape
k = math.log(m, diagLen)
if m != n or int(k) != k:
raise TypeError("The expression has shape {} so it cannot be "
"decomposed into subsystems of shape {}.".format(
glyphs.shape(self._shape), glyphs.shape((diagLen,)*2)))
return ((diagLen,)*2,)*int(k)
[docs] def partial_transpose(self, subsystems, dimensions=2):
r"""Return the expression with selected subsystems transposed.
If the expression can be written as
:math:`A_0 \otimes \cdots \otimes A_{n-1}` for matrices
:math:`A_0, \ldots, A_{n-1}` with shapes given in ``dimensions``, then
this returns :math:`B_0 \otimes \cdots \otimes B_{n-1}` with
:math:`B_i = A_i^T`, if ``i in subsystems`` (with :math:`i = -1` read as
:math:`n-1`), and :math:`B_i = A_i`, otherwise.
:param subsystems: A collection of or a single subystem number, indexed
from zero, corresponding to subsystems that shall be transposed.
The value :math:`-1` refers to the last subsystem.
:type subsystems: int or tuple or list
:param dimensions: Either an integer :math:`d` so that the subsystems
are assumed to be all of shape :math:`d \times d`, or a sequence of
subsystem shapes where an integer :math:`d` within the sequence is
read as :math:`d \times d`. In any case, the elementwise product
over all subsystem shapes must equal the expression's shape.
:type dimensions: int or tuple or list
:raises TypeError: If the subsystems do not match the expression.
:raises IndexError: If the subsystem selection is invalid.
:Example:
>>> from picos import Constant
>>> A = Constant("A", range(16), (4, 4))
>>> print(A) #doctest: +NORMALIZE_WHITESPACE
[ 0.00e+00 4.00e+00 8.00e+00 1.20e+01]
[ 1.00e+00 5.00e+00 9.00e+00 1.30e+01]
[ 2.00e+00 6.00e+00 1.00e+01 1.40e+01]
[ 3.00e+00 7.00e+00 1.10e+01 1.50e+01]
>>> A0 = A.partial_transpose(0); A0
<4×4 Real Constant: A.{[2×2]ᵀ⊗[2×2]}>
>>> print(A0) #doctest: +NORMALIZE_WHITESPACE
[ 0.00e+00 4.00e+00 2.00e+00 6.00e+00]
[ 1.00e+00 5.00e+00 3.00e+00 7.00e+00]
[ 8.00e+00 1.20e+01 1.00e+01 1.40e+01]
[ 9.00e+00 1.30e+01 1.10e+01 1.50e+01]
>>> A1 = A.partial_transpose(1); A1
<4×4 Real Constant: A.{[2×2]⊗[2×2]ᵀ}>
>>> print(A1) #doctest: +NORMALIZE_WHITESPACE
[ 0.00e+00 1.00e+00 8.00e+00 9.00e+00]
[ 4.00e+00 5.00e+00 1.20e+01 1.30e+01]
[ 2.00e+00 3.00e+00 1.00e+01 1.10e+01]
[ 6.00e+00 7.00e+00 1.40e+01 1.50e+01]
"""
m, n = self._shape
if isinstance(dimensions, int):
dimensions = self._square_equal_subsystem_dims(dimensions)
else:
dimensions = [
(d, d) if isinstance(d, int) else d for d in dimensions]
if reduce(
lambda x, y: (x[0]*y[0], x[1]*y[1]), dimensions) != self._shape:
raise TypeError("Subsystem dimensions do not match expression.")
if isinstance(subsystems, int):
subsystems = (subsystems,)
elif not subsystems:
return self
numSys = len(dimensions)
subsystems = set(numSys - 1 if sys == -1 else sys for sys in subsystems)
for sys in subsystems:
if not isinstance(sys, int):
raise IndexError("Subsystem indices must be integer, not {}."
.format(type(sys).__name__))
elif sys < 0:
raise IndexError("Subsystem indices must be nonnegative."
.format(sys))
elif sys >= numSys:
raise IndexError(
"Subsystem index {} out of range for {} systems total."
.format(sys, numSys))
# If all subsystems are transposed, this is regular transposition.
if len(subsystems) == numSys:
return self.T
# Prepare sparse K such that K·vec(A) = vec(partial_transpose(A)).
d = m * n
V = [1]*d
I = range(d)
J = cvxopt.matrix(I)
T = cvxopt.matrix(0, J.size)
obh, obw = 1, 1
sysStrings = None
# Apply transpositions starting with the rightmost Kronecker factor.
for sys in range(numSys - 1, -1, -1):
# Shape of current system.
p, q = dimensions[sys]
sysString = glyphs.matrix(glyphs.shape((p, q)))
# Height/width of "inner" blocks being moved, initially scalars.
ibh, ibw = obh, obw
# Heigh/width of "outer" blocks whose relative position is
# maintained but that are subject to transposition independently.
# In the last iteration this is the shape of the resulting matrix.
obh *= p
obw *= q
# Only transpose selected subsystems.
if sys not in subsystems:
sysStrings = glyphs.kron(sysString, sysStrings) \
if sysStrings else sysString
continue
else:
sysStrings = glyphs.kron(glyphs.transp(sysString), sysStrings) \
if sysStrings else glyphs.transp(sysString)
# Shape of outer blocks after transposition.
obhT, obwT = obw // ibw * ibh, obh // ibh * ibw
# Shape of full matrix after transposition.
mT, nT = m // obh * obhT, n // obw * obwT
for vi in I:
# Full matrix column and row.
c, r = divmod(vi, m)
# Outer block vertical index and row within outer block,
# outer block horizontal index and column within outer block.
obi, obr = divmod(r, obh)
obj, obc = divmod(c, obw)
# Inner block vertical index and row within inner block,
# inner block horizontal index and column within inner block.
ibi, ibr = divmod(obr, ibh)
ibj, ibc = divmod(obc, ibw)
# (1) ibi*ibw + ibc is column within the transposed outer block;
# adding obj*obwT yields the column in the transposed matrix.
# (2) ibj*ibh + ibr is row within the transposed outer block;
# adding obi*obhT yields the row in the transposed matrix.
# (3) tvi is index within the vectorized transposed matrix.
tvi = (obj*obwT + ibi*ibw + ibc)*mT \
+ (obi*obhT + ibj*ibh + ibr)
# Prepare the transposition.
T[tvi] = J[vi]
# Apply the transposition.
J, T = T, J
m, n, obh, obw = mT, nT, obhT, obwT
# Finalize the partial commutation matrix.
K = cvxopt.spmatrix(V, I, J, (d, d), self._typecode)
string = glyphs.ptransp_(self.string, sysStrings)
shape = (m, n)
coefs = {var: K * coef for var, coef in self._coefs.items()}
const = K * self._const
return self._basetype(string, shape, coefs, const)
@cached_property
def T0(self):
r"""Expression with the first :math:`2 \times 2` subsystem transposed.
Only available for a :math:`2^k \times 2^k` matrix with all subsystems
of shape :math:`2 \times 2`. Use :meth:`partial_transpose` otherwise.
"""
return self.partial_transpose(subsystems=0)
@cached_property
def T1(self):
r"""Expression with the second :math:`2 \times 2` subsystem transposed.
Only available for a :math:`2^k \times 2^k` matrix with all subsystems
of shape :math:`2 \times 2`. Use :meth:`partial_transpose` otherwise.
"""
return self.partial_transpose(subsystems=1)
@cached_property
def T2(self):
r"""Expression with the third :math:`2 \times 2` subsystem transposed.
Only available for a :math:`2^k \times 2^k` matrix with all subsystems
of shape :math:`2 \times 2`. Use :meth:`partial_transpose` otherwise.
"""
return self.partial_transpose(subsystems=2)
@cached_property
def T3(self):
r"""Expression with the fourth :math:`2 \times 2` subsystem transposed.
Only available for a :math:`2^k \times 2^k` matrix with all subsystems
of shape :math:`2 \times 2`. Use :meth:`partial_transpose` otherwise.
"""
return self.partial_transpose(subsystems=3)
@cached_property
def Tl(self):
r"""Expression with the last :math:`2 \times 2` subsystem transposed.
Only available for a :math:`2^k \times 2^k` matrix with all subsystems
of shape :math:`2 \times 2`. Use :meth:`partial_transpose` otherwise.
"""
return self.partial_transpose(subsystems=-1)
@staticmethod
def _reindex_F(indices, source, destination):
"""Convert indices between different tensor shapes in Fortran-order."""
new = []
offset = 0
factor = 1
for index, base in zip(indices, source):
offset += factor*index
factor *= base
for base in destination:
offset, remainder = divmod(offset, base)
new.append(remainder)
return tuple(new)
@staticmethod
def _reindex_C(indices, source, destination):
"""Convert indices between different tensor shapes in C-order."""
new = []
offset = 0
factor = 1
for index, base in zip(reversed(indices), reversed(source)):
offset += factor*index
factor *= base
for base in reversed(destination):
offset, remainder = divmod(offset, base)
new.insert(0, remainder)
return tuple(new)
[docs] def reshuffled(self, permutation="ikjl", dimensions=None, order="C"):
"""Return the reshuffled or realigned expression.
This operation works directly on matrices. However, it is equivalent to
the following sequence of operations:
1. The matrix is reshaped to a tensor with the given ``dimensions`` and
according to ``order``.
2. The tensor's axes are permuted according to ``permutation``.
3. The tensor is reshaped back to the shape of the original matrix
according to ``order``.
For comparison, the following function applies the same operation to a
2D NumPy :class:`~numpy:numpy.ndarray`:
.. code::
def reshuffle_numpy(matrix, permutation, dimensions, order):
P = "{} -> {}".format("".join(sorted(permutation)), permutation)
reshuffled = numpy.reshape(matrix, dimensions, order)
reshuffled = numpy.einsum(P, reshuffled)
return numpy.reshape(reshuffled, matrix.shape, order)
:param permutation:
A sequence of comparable elements with length equal to the number of
tensor dimensions. The sequence is compared to its ordered version
and the resulting permutation pattern is used to permute the tensor
indices. For instance, the string ``"ikjl"`` is compared to its
sorted version ``"ijkl"`` and denotes that the second and third axis
should be swapped.
:type permutation:
str or tuple or list
:param dimensions:
If this is an integer sequence, then it defines the dimensions of
the tensor. If this is :obj:`None`, then the tensor is assumed to be
hypercubic and the number of dimensions is inferred from the
``permutation`` argument.
:type dimensions:
None or tuple or list
:param str order:
The indexing order to use for the virtual reshaping. Must be either
``"F"`` for Fortran-order (generalization of column-major) or
``"C"`` for C-order (generalization of row-major). Note that PICOS
usually reshapes in Fortran-order while NumPy defaults to C-order.
:Example:
>>> from picos import Constant
>>> A = Constant("A", range(16), (4, 4))
>>> print(A) #doctest: +NORMALIZE_WHITESPACE
[ 0.00e+00 4.00e+00 8.00e+00 1.20e+01]
[ 1.00e+00 5.00e+00 9.00e+00 1.30e+01]
[ 2.00e+00 6.00e+00 1.00e+01 1.40e+01]
[ 3.00e+00 7.00e+00 1.10e+01 1.50e+01]
>>> R = A.reshuffled(); R
<4×4 Real Constant: shuffled(A,ikjl,C)>
>>> print(R) #doctest: +NORMALIZE_WHITESPACE
[ 0.00e+00 4.00e+00 1.00e+00 5.00e+00]
[ 8.00e+00 1.20e+01 9.00e+00 1.30e+01]
[ 2.00e+00 6.00e+00 3.00e+00 7.00e+00]
[ 1.00e+01 1.40e+01 1.10e+01 1.50e+01]
>>> A.reshuffled("ji").equals(A.T) # Regular transposition.
True
>>> A.reshuffled("3214").equals(A.T0) # Partial transposition (1).
True
>>> A.reshuffled("1432").equals(A.T1) # Partial transposition (2).
True
"""
m, n = self._shape
mn = m*n
# Load the permutation.
ordered = sorted(permutation)
P = dict(enumerate(ordered.index(element) for element in permutation))
if len(set(P.values())) < len(P):
raise ValueError("The sequence defining the permutation appears to "
"contain duplicate elements.")
assert not set(P.keys()).symmetric_difference(set(P.values()))
numDims = len(P)
# Load the dimensions.
guessDimensions = dimensions is None
if guessDimensions:
dimensions = (int(mn**(1.0 / numDims)),)*numDims
else:
if len(dimensions) != numDims:
raise ValueError("The number of indices does not match the "
"number of dimensions.")
if reduce(int.__mul__, dimensions, 1) != mn:
raise TypeError("The {} matrix {} cannot be reshaped to a {} "
"tensor.".format(glyphs.shape(self.shape), self.string,
"hypercubic order {}".format(numDims) if guessDimensions
else glyphs.size("", "").join(str(d) for d in dimensions)))
# Load the indexing order.
if order not in "FC":
raise ValueError("Order must be given as 'F' or 'C'.")
reindex = self._reindex_F if order == "F" else self._reindex_C
# Nothing to do for the neutral permutation.
if all(key == val for key, val in P.items()):
return self
# Create a sparse mtrix R such that R·vec(A) = vec(reshuffle(A)).
V, I, J = [1]*mn, [], range(mn)
for i in range(mn):
(k, j) = divmod(i, m) # (j, k) are column-major matrix indices.
indices = reindex((j, k), (m, n), dimensions)
newIndices = tuple(indices[P[d]] for d in range(numDims))
newDimensions = tuple(dimensions[P[d]] for d in range(numDims))
(j, k) = reindex(newIndices, newDimensions, (m, n))
I.append(k*m + j)
R = cvxopt.spmatrix(V, I, J, (mn, mn), self._typecode)
# Create the string.
strArgs = [self.string, str(permutation).replace(" ", ""), order]
if not guessDimensions:
strArgs.insert(2, str(dimensions).replace(" ", ""))
string = glyphs.shuffled(",".join(strArgs))
# Finalize the new expression.
shape = (m, n)
coefs = {var: R * coef for var, coef in self._coefs.items()}
const = R * self._const
return self._basetype(string, shape, coefs, const)
@cached_property
def sum(self):
"""Sum over all scalar elements of the expression."""
# NOTE: glyphs.clever_dotp detects this case and uses the sum glyph.
# NOTE: 1 on the right hand side in case self is complex.
return (self | 1)
@cached_property
def tr(self):
"""Trace of a square expression."""
if not self.square:
raise TypeError("Cannot compute the trace of non-square {}."
.format(self.string))
# NOTE: glyphs.clever_dotp detects this case and uses the trace glyph.
# NOTE: "I" on the right hand side in case self is complex.
return (self | "I")
[docs] def partial_trace(self, subsystems, dimensions=2):
r"""Return the partial trace over selected subsystems.
If the expression can be written as
:math:`A_0 \otimes \cdots \otimes A_{n-1}` for matrices
:math:`A_0, \ldots, A_{n-1}` with shapes given in ``dimensions``, then
this returns :math:`B_0 \otimes \cdots \otimes B_{n-1}` with
:math:`B_i = \operatorname{tr}(A_i)`, if ``i in subsystems`` (with
:math:`i = -1` read as :math:`n-1`), and :math:`B_i = A_i`, otherwise.
:param subsystems: A collection of or a single subystem number, indexed
from zero, corresponding to subsystems that shall be traced over.
The value :math:`-1` refers to the last subsystem.
:type subsystems: int or tuple or list
:param dimensions: Either an integer :math:`d` so that the subsystems
are assumed to be all of shape :math:`d \times d`, or a sequence of
subsystem shapes where an integer :math:`d` within the sequence is
read as :math:`d \times d`. In any case, the elementwise product
over all subsystem shapes must equal the expression's shape.
:type dimensions: int or tuple or list
:raises TypeError: If the subsystems do not match the expression or if
a non-square subsystem is to be traced over.
:raises IndexError: If the subsystem selection is invalid in any other
way.
:Example:
>>> from picos import Constant
>>> A = Constant("A", range(16), (4, 4))
>>> print(A) #doctest: +NORMALIZE_WHITESPACE
[ 0.00e+00 4.00e+00 8.00e+00 1.20e+01]
[ 1.00e+00 5.00e+00 9.00e+00 1.30e+01]
[ 2.00e+00 6.00e+00 1.00e+01 1.40e+01]
[ 3.00e+00 7.00e+00 1.10e+01 1.50e+01]
>>> A0 = A.partial_trace(0); A0
<2×2 Real Constant: A.{tr([2×2])⊗[2×2]}>
>>> print(A0) #doctest: +NORMALIZE_WHITESPACE
[ 1.00e+01 1.80e+01]
[ 1.20e+01 2.00e+01]
>>> A1 = A.partial_trace(1); A1
<2×2 Real Constant: A.{[2×2]⊗tr([2×2])}>
>>> print(A1) #doctest: +NORMALIZE_WHITESPACE
[ 5.00e+00 2.10e+01]
[ 9.00e+00 2.50e+01]
"""
# Shape of the original matrix.
m, n = self._shape
if isinstance(dimensions, int):
dimensions = self._square_equal_subsystem_dims(dimensions)
else:
dimensions = [
(d, d) if isinstance(d, int) else d for d in dimensions]
if reduce(
lambda x, y: (x[0]*y[0], x[1]*y[1]), dimensions) != self._shape:
raise TypeError("Subsystem dimensions do not match expression.")
if isinstance(subsystems, int):
subsystems = (subsystems,)
elif not subsystems:
return self
numSys = len(dimensions)
subsystems = set(numSys - 1 if sys == -1 else sys for sys in subsystems)
for sys in subsystems:
if not isinstance(sys, int):
raise IndexError("Subsystem indices must be integer, not {}."
.format(type(sys).__name__))
elif sys < 0:
raise IndexError("Subsystem indices must be nonnegative."
.format(sys))
elif sys >= numSys:
raise IndexError(
"Subsystem index {} out of range for {} systems total."
.format(sys, numSys))
elif dimensions[sys][0] != dimensions[sys][1]:
raise TypeError(
"Subsystem index {} refers to a non-square subsystem that "
"cannot be traced over.".format(sys))
# If all subsystems are traced over, this is the regular trace.
if len(subsystems) == numSys:
return self.tr
# Prepare sparse T such that T·vec(A) = vec(partial_trace(A)).
T = []
# Compute factors of T, one for each subsystem being traced over.
ibh, ibw = m, n
sysStrings = None
for sys in range(numSys):
# Shape of current system.
p, q = dimensions[sys]
sysString = glyphs.matrix(glyphs.shape((p, q)))
# Heigh/width of "outer" blocks whose relative position is
# maintained but that are reduced independently to the size of an
# inner block if the current system is to be traced over.
obh, obw = ibh, ibw
# Height/width of "inner" blocks that are summed if they are
# main-diagonal blocks of an outer block.
ibh, ibw = obh // p, obw // q
# Only trace over selected subsystems.
if sys not in subsystems:
sysStrings = glyphs.kron(sysStrings, sysString) \
if sysStrings else sysString
continue
else:
sysStrings = glyphs.kron(sysStrings, glyphs.trace(sysString)) \
if sysStrings else glyphs.trace(sysString)
# Shape of new matrix.
assert p == q
mN, nN = m // p, n // p
# Prepare one factor of T.
oldLen = m * n
newLen = mN * nN
V, I, J = [1]*(newLen*p), [], []
shape = (newLen, oldLen)
# Determine the summands that make up each entry of the new matrix.
for viN in range(newLen):
# A column/row pair that represents a scalar in the new matrix
# and a number of p scalars within different on-diagonal inner
# blocks but within the same outer block in the old matrix.
cN, rN = divmod(viN, mN)
# Index pair (obi, obj) for outer block in question, row/column
# pair (ibr, ibc) identifies the scalar within each inner block.
obi, ibr = divmod(rN, ibh)
obj, ibc = divmod(cN, ibw)
# Collect summands for the current entry of the new matrix; one
# scalar per on-diagonal inner block.
for k in range(p):
rO = obi*obh + k*ibh + ibr
cO = obj*obw + k*ibw + ibc
I.append(viN)
J.append(cO*m + rO)
# Store the operator that performs the current subsystem trace.
T.insert(0, cvxopt.spmatrix(V, I, J, shape, self._typecode))
# Make next iteration work on the new matrix.
m, n = mN, nN
# Multiply out the linear partial trace operator T.
# TODO: Fast matrix multiplication dynamic program useful here?
T = reduce(lambda A, B: A*B, T)
string = glyphs.ptrace_(self.string, sysStrings)
shape = (m, n)
coefs = {var: T * coef for var, coef in self._coefs.items()}
const = T * self._const
return self._basetype(string, shape, coefs, const)
@cached_property
def tr0(self):
r"""Expression with the first :math:`2 \times 2` subsystem traced out.
Only available for a :math:`2^k \times 2^k` matrix with all subsystems
of shape :math:`2 \times 2`. Use :meth:`partial_trace` otherwise.
"""
return self.partial_trace(subsystems=0)
@cached_property
def tr1(self):
r"""Expression with the second :math:`2 \times 2` subsystem traced out.
Only available for a :math:`2^k \times 2^k` matrix with all subsystems
of shape :math:`2 \times 2`. Use :meth:`partial_trace` otherwise.
"""
return self.partial_trace(subsystems=1)
@cached_property
def tr2(self):
r"""Expression with the third :math:`2 \times 2` subsystem traced out.
Only available for a :math:`2^k \times 2^k` matrix with all subsystems
of shape :math:`2 \times 2`. Use :meth:`partial_trace` otherwise.
"""
return self.partial_trace(subsystems=2)
@cached_property
def tr3(self):
r"""Expression with the fourth :math:`2 \times 2` subsystem traced out.
Only available for a :math:`2^k \times 2^k` matrix with all subsystems
of shape :math:`2 \times 2`. Use :meth:`partial_trace` otherwise.
"""
return self.partial_trace(subsystems=3)
@cached_property
def trl(self):
r"""Expression with the last :math:`2 \times 2` subsystem traced out.
Only available for a :math:`2^k \times 2^k` matrix with all subsystems
of shape :math:`2 \times 2`. Use :meth:`partial_trace` otherwise.
"""
return self.partial_trace(subsystems=-1)
@cached_property
def vec(self):
"""Column-major vectorization of the expression as a column vector.
.. note::
Given an expression ``A``, ``A.vec`` and ``A[:]`` produce the same
result (up to its string description) but ``A.vec`` is faster and
its result is cached.
:Example:
>>> from picos import Constant
>>> A = Constant("A", [[1, 2], [3, 4]])
>>> A.vec.equals(A[:])
True
>>> A[:] is A[:]
False
>>> A.vec is A.vec
True
"""
if self._shape[1] == 1:
return self
else:
return self._basetype(glyphs.vec(self.string), (len(self), 1),
self._coefs, self._const)
[docs] def dupvec(self, n):
"""Return a (repeated) column-major vectorization of the expression.
:param int n: Number of times to duplicate the vectorization.
:returns: A column vector.
:Example:
>>> from picos import Constant
>>> A = Constant("A", [[1, 2], [3, 4]])
>>> A.dupvec(1) is A.vec
True
>>> A.dupvec(3).equals(A.vec // A.vec // A.vec)
True
"""
if not isinstance(n, int):
raise TypeError("Number of copies must be integer.")
if n < 1:
raise ValueError("Number of copies must be positive.")
if n == 1:
return self.vec
else:
string = glyphs.vec(glyphs.comma(self.string, n))
shape = (len(self)*n, 1)
coefs = {var: cvxopt_vcat([coef]*n)
for var, coef in self._coefs.items()}
const = cvxopt_vcat([self._const]*n)
return self._basetype(string, shape, coefs, const)
@cached_property
def trilvec(self):
r"""Column-major vectorization of the lower triangular part.
:returns:
A column vector of all elements :math:`A_{ij}` that satisfy
:math:`i \geq j`.
.. note::
If you want a row-major vectorization instead, write ``A.T.triuvec``
instead of ``A.trilvec``.
:Example:
>>> from picos import Constant
>>> A = Constant("A", [[1, 2], [3, 4], [5, 6]])
>>> print(A)
[ 1.00e+00 2.00e+00]
[ 3.00e+00 4.00e+00]
[ 5.00e+00 6.00e+00]
>>> print(A.trilvec)
[ 1.00e+00]
[ 3.00e+00]
[ 5.00e+00]
[ 4.00e+00]
[ 6.00e+00]
"""
m, n = self._shape
if n == 1: # Column vector or scalar.
return self
elif m == 1: # Row vector.
return self[0]
# Build a transformation D such that D·vec(A) = trilvec(A).
rows = [j*m + i for j in range(n) for i in range(m) if i >= j]
d = len(rows)
D = cvxopt.spmatrix([1]*d, range(d), rows, (d, len(self)))
string = glyphs.trilvec(self.string)
shape = (d, 1)
coefs = {var: D*coef for var, coef in self._coefs.items()}
const = D*self._const
return self._basetype(string, shape, coefs, const)
@cached_property
def triuvec(self):
r"""Column-major vectorization of the upper triangular part.
:returns:
A column vector of all elements :math:`A_{ij}` that satisfy
:math:`i \leq j`.
.. note::
If you want a row-major vectorization instead, write ``A.T.trilvec``
instead of ``A.triuvec``.
:Example:
>>> from picos import Constant
>>> A = Constant("A", [[1, 2, 3], [4, 5, 6]])
>>> print(A)
[ 1.00e+00 2.00e+00 3.00e+00]
[ 4.00e+00 5.00e+00 6.00e+00]
>>> print(A.triuvec)
[ 1.00e+00]
[ 2.00e+00]
[ 5.00e+00]
[ 3.00e+00]
[ 6.00e+00]
"""
m, n = self._shape
if m == 1: # Row vector or scalar.
return self
elif n == 1: # Column vector.
return self[0]
# Build a transformation D such that D·vec(A) = triuvec(A).
rows = [j*m + i for j in range(n) for i in range(m) if i <= j]
d = len(rows)
D = cvxopt.spmatrix([1]*d, range(d), rows, (d, len(self)))
string = glyphs.triuvec(self.string)
shape = (d, 1)
coefs = {var: D*coef for var, coef in self._coefs.items()}
const = D*self._const
return self._basetype(string, shape, coefs, const)
[docs] def dupdiag(self, n):
"""Return a matrix with the (repeated) expression on the diagonal.
Vectorization is performed in column-major order.
:param int n: Number of times to duplicate the vectorization.
"""
# Vectorize and duplicate the expression.
vec = self.dupvec(n)
d = len(vec)
# Build a transformation D such that D·vec(A) = vec(diag(vec(A))).
ones = [1]*d
D = cvxopt.spdiag(ones)[:]
D = cvxopt.spmatrix(ones, D.I, range(d), (D.size[0], d))
string = glyphs.diag(vec.string)
shape = (d, d)
coefs = {var: D*coef for var, coef in vec._coefs.items()}
const = D*vec._const
return self._basetype(string, shape, coefs, const)
@cached_property
def diag(self):
"""Diagonal matrix with the expression on the main diagonal.
Vectorization is performed in column-major order.
"""
return self.dupdiag(1)
@cached_property
def maindiag(self):
"""The main diagonal of the expression as a column vector."""
if 1 in self._shape:
return self[0]
# Build a transformation D such that D·vec(A) = diag(A).
step = self._shape[0] + 1
rows = [i*step for i in range(min(self._shape))]
d = len(rows)
D = cvxopt.spmatrix([1]*d, range(d), rows, (d, len(self)))
string = glyphs.maindiag(self.string)
shape = (d, 1)
coefs = {var: D*coef for var, coef in self._coefs.items()}
const = D*self._const
return self._basetype(string, shape, coefs, const)
# --------------------------------------------------------------------------
# Constraint-creating operators, and _predict.
# --------------------------------------------------------------------------
@classmethod
def _predict(cls, subtype, relation, other):
assert isinstance(subtype, cls.Subtype)
from .set import Set
if relation == operator.__eq__:
if issubclass(other.clstype, ComplexAffineExpression):
return ComplexAffineConstraint.make_type(dim=subtype.dim)
elif relation == operator.__lshift__:
if issubclass(other.clstype, ComplexAffineExpression):
return ComplexLMIConstraint.make_type(int(subtype.dim**0.5))
elif issubclass(other.clstype, Set):
return other >> ExpressionType(cls, subtype)
elif relation == operator.__rshift__:
if issubclass(other.clstype, ComplexAffineExpression):
return ComplexLMIConstraint.make_type(int(subtype.dim**0.5))
return NotImplemented
@convert_operands(sameShape=True)
@validate_prediction
@refine_operands()
def __eq__(self, other):
if isinstance(other, ComplexAffineExpression):
return ComplexAffineConstraint(self, other)
else:
return NotImplemented
# Since we define __eq__, __hash__ is not inherited. Do this manually.
__hash__ = Expression.__hash__
@convert_operands(sameShape=True)
@validate_prediction
@refine_operands()
def __lshift__(self, other):
from .set import Set
if isinstance(other, ComplexAffineExpression):
return ComplexLMIConstraint(self, Constraint.LE, other)
elif isinstance(other, Set):
# NOTE: We need to swap explicitly because Python does not consider
# A << B the same as B >> A.
return other.__rshift__(self)
else:
return NotImplemented
@convert_operands(sameShape=True)
@validate_prediction
@refine_operands()
def __rshift__(self, other):
if isinstance(other, ComplexAffineExpression):
return ComplexLMIConstraint(self, Constraint.GE, other)
else:
return NotImplemented
# --------------------------------------------------------------------------
# Advanced interface for PICOS-internal use.
# --------------------------------------------------------------------------
[docs] def sparse_rows(
self, varOffsetMap, lowerTriangle=False, upperTriangle=False,
indexFunction=None):
"""Return a sparse list representation of the expression.
The method is intended for internal use: It simplifies passing affine
constraints to solvers that support only scalar constraints. The idea is
to pose the constraint as a single (multidimensional) affine expression
bounded by zero, and use the coefficients and the constant term of this
expression to fill the solver's constraint matrix (with columns
representing scalar variables and rows representing scalar constraints).
:param dict varOffsetMap: Maps variables to column offsets.
:param bool lowerTriangle: Whether to return only the lower triangular
part of the expression.
:param bool upperTriangle: Whether to return only the upper triangular
part of the expression.
:param indexFunction: Instead of adding the local variable index to the
value returned by varOffsetMap, use the return value of this
function, that takes as argument the variable and its local index,
as the "column index", which need not be an integer. When this
parameter is passed, the parameter varOffsetMap is ignored.
:returns: A list of triples (J, V, c) where J contains column indices
(representing scalar variables), V contains coefficients for each
column index and c is a constant term.
"""
if lowerTriangle and upperTriangle:
lowerTriangle = False
upperTriangle = False
rows = []
numRows = len(self)
m = self.size[0]
for row in range(numRows):
i, j = row % m, row // m
if lowerTriangle and i < j:
rows.append(None)
elif upperTriangle and j < i:
rows.append(None)
else:
rows.append([[], [], 0.0])
for var, coef in self._coefs.items():
V, I, J, _ = sparse_quadruple(coef)
for localCoef, localConIndex, localVarIndex in zip(V, I, J):
row = rows[localConIndex]
if not row:
continue
# TODO: Use a single parameter for both types.
if indexFunction:
row[0].append(indexFunction(var, localVarIndex))
else:
row[0].append(varOffsetMap[var] + localVarIndex)
row[1].append(localCoef)
for localConIndex in range(numRows):
row = rows[localConIndex]
if not row:
continue
row[2] = self._const[localConIndex]
return [row for row in rows if row]
# --------------------------------------------------------------------------
# Backwards compatibility methods.
# --------------------------------------------------------------------------
[docs] @classmethod
@deprecated("2.0", useInstead="from_constant")
def fromScalar(cls, scalar):
"""Create a :class:`ComplexAffineExpression` from a numeric scalar."""
return cls.from_constant(scalar, (1, 1))
[docs] @classmethod
@deprecated("2.0", useInstead="from_constant")
def fromMatrix(cls, matrix, size=None):
"""Create a :class:`ComplexAffineExpression` from a numeric matrix."""
return cls.from_constant(matrix, size)
[docs] @deprecated("2.0", useInstead="object.__xor__")
def hadamard(self, fact):
"""Denote the elementwise (or Hadamard) product."""
return self.__xor__(fact)
[docs] @deprecated("2.0", useInstead="constant")
def isconstant(self):
"""Whether the expression involves no variables."""
return not self._coefs
[docs] @deprecated("2.0", useInstead="equals")
def same_as(self, other):
"""Check mathematical equality with another affine expression."""
return self.equals(other)
[docs] @deprecated("2.0", useInstead="T")
def transpose(self):
"""Return the matrix transpose."""
return self.T
@cached_property
@deprecated("2.0", useInstead="partial_transpose", decoratorLevel=1)
def Tx(self):
"""Auto-detect few subsystems of same shape and transpose the last."""
m, n = self._shape
dims = None
for k in range(2, int(math.log(min(m, n), 2)) + 1):
p, q = int(round(m**(1.0/k))), int(round(n**(1.0/k)))
if m == p**k and n == q**k:
dims = ((p, q),)*k
break
if dims:
return self.partial_transpose(subsystems=-1, dimensions=dims)
else:
raise RuntimeError("Failed to auto-detect subsystem dimensions for "
"partial transposition: Only supported for {} matrices, {}."
.format(glyphs.shape(
(glyphs.power("m", "k"), glyphs.power("n", "k"))),
glyphs.ge("k", 2)))
[docs] @deprecated("2.0", useInstead="conj")
def conjugate(self):
"""Return the complex conjugate."""
return self.conj
[docs] @deprecated("2.0", useInstead="H")
def Htranspose(self):
"""Return the conjugate (or Hermitian) transpose."""
return self.H
[docs] @deprecated("2.0", reason="PICOS expressions are now immutable.")
def copy(self):
"""Return a deep copy of the expression."""
from copy import copy as cp
return self._basetype(
cp(self._symbStr),
self._shape,
{var: cp(coef) for var, coef in self._coefs.items()},
cp(self._const))
[docs] @deprecated("2.0", reason="PICOS expressions are now immutable.")
def soft_copy(self):
"""Return a shallow copy of the expression."""
return self._basetype(
self._symbStr, self._shape, self._coefs, self._const)
[docs]class AffineExpression(ComplexAffineExpression):
"""A multidimensional real affine expression."""
# --------------------------------------------------------------------------
# Method overridings.
# --------------------------------------------------------------------------
@classmethod
def _get_basetype(cls):
return AffineExpression
@classmethod
def _get_typecode(cls):
return "d"
@classmethod
def _get_type_string_base(cls):
return "Real {}"
def _get_refined(self):
return self
@convert_operands(sameShape=True, allowNone=True)
def _set_value(self, value):
if value is None:
for var in self._coefs:
var.value = None
return
if not isinstance(value, AffineExpression) or not value.constant:
raise TypeError("Cannot set the value of {} to {}: Not real or not "
"a constant.".format(repr(self), repr(value)))
if self.constant:
raise TypeError("Cannot set the value on a constant expression.")
y = cvx2np(value._const)
A = []
for var, coef in self._coefs.items():
A.append(cvx2np(coef))
assert A
A = numpy.hstack(A)
b = y - cvx2np(self._const)
try:
solution, residual, _, _ = numpy.linalg.lstsq(A, b, rcond=None)
except numpy.linalg.LinAlgError as error:
raise RuntimeError("Setting a value on {} failed: {}."
.format(self.string, error))
if not numpy.allclose(residual, 0):
raise ValueError("Setting a value on {} failed: No exact solution "
"to the associated linear system found.".format(self.string))
offset = 0
for var in self._coefs:
var.internal_value = solution[offset:offset+var.dim]
offset += var.dim
@property
def isreal(self):
"""Always true for :class:`AffineExpression` instances.""" # noqa
return True
@property
def real(self):
"""The :class:`AffineExpression` as is.""" # noqa
return self
@cached_property
def imag(self):
"""A zero of same shape as the :class:`AffineExpression`.""" # noqa
return self._basetype.zero(self._shape)
@property
def conj(self):
"""The :class:`AffineExpression` as is.""" # noqa
return self
@property
def H(self):
"""The regular transpose of the :class:`AffineExpression`.""" # noqa
return self.T
# --------------------------------------------------------------------------
# Additional methods for real affine expressions only.
# --------------------------------------------------------------------------
@cached_property
def exp(self):
"""The exponential function applied to the expression.""" # noqa
from . import SumExponentials
return SumExponentials(self)
@cached_property
def log(self):
"""The Logarithm of the expression.""" # noqa
from . import Logarithm
return Logarithm(self)
# --------------------------------------------------------------------------
# Constraint-creating operators, and _predict.
# --------------------------------------------------------------------------
@classmethod
def _predict(cls, subtype, relation, other):
assert isinstance(subtype, cls.Subtype)
from .set import Set
if relation in (operator.__eq__, operator.__le__, operator.__ge__):
if issubclass(other.clstype, AffineExpression):
return AffineConstraint.make_type(
dim=subtype.dim, eq=(relation is operator.__eq__))
elif relation == operator.__lshift__:
if issubclass(other.clstype, AffineExpression):
return LMIConstraint.make_type(int(subtype.dim**0.5))
elif issubclass(other.clstype, ComplexAffineExpression):
return ComplexLMIConstraint.make_type(int(subtype.dim**0.5))
elif issubclass(other.clstype, Set):
return other >> ExpressionType(cls, subtype)
elif relation == operator.__rshift__:
if issubclass(other.clstype, AffineExpression):
return LMIConstraint.make_type(int(subtype.dim**0.5))
elif issubclass(other.clstype, ComplexAffineExpression):
return ComplexLMIConstraint.make_type(int(subtype.dim**0.5))
return NotImplemented
@convert_operands(sameShape=True)
@validate_prediction
@refine_operands()
def __le__(self, other):
if isinstance(other, AffineExpression):
return AffineConstraint(self, Constraint.LE, other)
else:
return NotImplemented
@convert_operands(sameShape=True)
@validate_prediction
@refine_operands()
def __ge__(self, other):
if isinstance(other, AffineExpression):
return AffineConstraint(self, Constraint.GE, other)
else:
return NotImplemented
@convert_operands(sameShape=True)
@validate_prediction
@refine_operands()
def __eq__(self, other):
if isinstance(other, AffineExpression):
return AffineConstraint(self, Constraint.EQ, other)
else:
return NotImplemented
# Since we define __eq__, __hash__ is not inherited. Do this manually.
__hash__ = Expression.__hash__
@convert_operands(sameShape=True)
@validate_prediction
@refine_operands()
def __lshift__(self, other):
from .set import Set
if isinstance(other, AffineExpression):
return LMIConstraint(self, Constraint.LE, other)
elif isinstance(other, ComplexAffineExpression):
return ComplexLMIConstraint(self, Constraint.LE, other)
elif isinstance(other, Set):
# NOTE: We need to swap explicitly because Python does not consider
# A << B the same as B >> A.
return other.__rshift__(self)
else:
return NotImplemented
@convert_operands(sameShape=True)
@validate_prediction
@refine_operands()
def __rshift__(self, other):
if isinstance(other, AffineExpression):
return LMIConstraint(self, Constraint.GE, other)
elif isinstance(other, ComplexAffineExpression):
return ComplexLMIConstraint(self, Constraint.GE, other)
else:
return NotImplemented
[docs]def Constant(name_or_value, value=None, shape=None):
"""Create a constant PICOS expression.
Loads the given numeric value as a constant
:class:`~picos.expressions.ComplexAffineExpression` or
:class:`~picos.expressions.AffineExpression`, depending on the value.
Optionally, the value is broadcasted or reshaped according to the shape
argument.
:param str name_or_value: Symbolic string description of the constant. If
:obj:`None` or the empty string, a string will be generated. If this is
the only positional parameter (i.e.``value`` is not given), then this
position is used as the value argument instead!
:param value: The numeric constant to load.
See :func:`~.data.load_data` for supported data formats and broadcasting and
reshaping rules.
:Example:
>>> from picos import Constant
>>> Constant(1)
<1×1 Real Constant: 1>
>>> Constant(1, shape=(2, 2))
<2×2 Real Constant: [1]>
>>> Constant("one", 1)
<1×1 Real Constant: one>
>>> Constant("J", 1, (2, 2))
<2×2 Real Constant: J>
"""
if value is None:
value = name_or_value
name = None
else:
name = name_or_value
value, valStr = load_data(value, shape)
if value.typecode == "z":
cls = ComplexAffineExpression
else:
cls = AffineExpression
return cls(name if name else valStr, value.size, {}, value)
# --------------------------------------
__all__ = api_end(_API_START, globals())