Coverage for picos/expressions/exp_affine.py: 85.71%
315 statements
« prev ^ index » next coverage.py v6.5.0, created at 2023-03-26 07:46 +0000
1# ------------------------------------------------------------------------------
2# Copyright (C) 2019-2022 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 affine expression types."""
22import operator
23from collections import namedtuple
25import cvxopt
26import numpy
28from .. import glyphs
29from ..apidoc import api_end, api_start
30from ..caching import cached_property, cached_unary_operator
31from ..constraints import (AffineConstraint, ComplexAffineConstraint,
32 ComplexLMIConstraint, Constraint, LMIConstraint)
33from .data import convert_operands, cvx2np, cvx2csc, load_data
34from .exp_biaffine import BiaffineExpression
35from .expression import (Expression, ExpressionType, refine_operands,
36 validate_prediction)
38_API_START = api_start(globals())
39# -------------------------------
42class ComplexAffineExpression(BiaffineExpression):
43 """A multidimensional (complex) affine expression.
45 Base class for the real :class:`AffineExpression`.
46 """
48 # --------------------------------------------------------------------------
49 # Abstract method implementations for Expression.
50 # --------------------------------------------------------------------------
52 Subtype = namedtuple("Subtype", ("shape", "constant", "nonneg"))
53 Subtype.dim = property(lambda self: self.shape[0] * self.shape[1])
55 def _get_subtype(self):
56 """Implement :meth:`~.expression.Expression._get_subtype`."""
57 nonneg = self.constant and self.isreal \
58 and all(x >= 0 for x in self.value_as_matrix)
60 return self.Subtype(self._shape, self.constant, nonneg)
62 # --------------------------------------------------------------------------
63 # Method overridings for Expression.
64 # --------------------------------------------------------------------------
66 @cached_unary_operator
67 def _get_refined(self):
68 """Implement :meth:`~.expression.Expression._get_refined`."""
69 if self.isreal:
70 return AffineExpression(self._symbStr, self._shape, self._coefs)
71 else:
72 return self
74 @convert_operands(sameShape=True, allowNone=True)
75 def _set_value(self, value):
76 """Override :meth:`~.expression.Expression._set_value`."""
77 if value is None:
78 for var in self._linear_coefs:
79 var.value = None
80 return
82 # Since all variables are real-valued, prevent NumPy from finding
83 # complex-valued solutions that do not actually work.
84 (self.real // self.imag).renamed(self.string).value \
85 = (value.real // value.imag)
87 # --------------------------------------------------------------------------
88 # Abstract method implementations for BiaffineExpression.
89 # --------------------------------------------------------------------------
91 @classmethod
92 def _get_bilinear_terms_allowed(cls):
93 """Implement for :class:`~.exp_biaffine.BiaffineExpression`."""
94 return False
96 @classmethod
97 def _get_parameters_allowed(cls):
98 """Implement for :class:`~.exp_biaffine.BiaffineExpression`."""
99 return False
101 @classmethod
102 def _get_basetype(cls):
103 """Implement :meth:`~.exp_biaffine.BiaffineExpression._get_basetype`."""
104 return ComplexAffineExpression
106 @classmethod
107 def _get_typecode(cls):
108 """Implement :meth:`~.exp_biaffine.BiaffineExpression._get_typecode`."""
109 return "z"
111 # --------------------------------------------------------------------------
112 # Method overridings for BiaffineExpression: Binary operators.
113 # --------------------------------------------------------------------------
115 @convert_operands(sameShape=True)
116 @refine_operands(stop_at_affine=True)
117 def __or__(self, other):
118 from .exp_quadratic import QuadraticExpression
119 from .exp_sqnorm import SquaredNorm
121 if isinstance(other, ComplexAffineExpression) \
122 and not self.constant and not other.constant:
123 # Create a squared norm if possible.
124 # NOTE: Must not check self.equals(other) here; see SquaredNorm.
125 # TODO: Consider creating a helper function for __or__ that always
126 # returns a QuadraticExpression instead of a SquaredNorm to be
127 # used within SquaredNorm. Then equals would be possible here.
128 if self is other:
129 return SquaredNorm(self)
131 string = glyphs.clever_dotp(
132 self.string, other.string, other.complex, self.scalar)
134 # Handle the complex case: Conjugate the right hand side.
135 other = other.conj
137 Cs, Co = self._constant_coef, other._constant_coef
139 # Compute the affine part of the product.
140 affString = glyphs.affpart(string)
141 affCoefs = {(): Cs.T * Co}
142 for var in self.variables.union(other.variables):
143 if var not in other._linear_coefs:
144 affCoefs[var] = Co.T * self._linear_coefs[var]
145 elif var not in self._linear_coefs:
146 affCoefs[var] = Cs.T * other._linear_coefs[var]
147 else:
148 affCoefs[var] = Co.T * self._linear_coefs[var] + \
149 Cs.T * other._linear_coefs[var]
150 affPart = self._common_basetype(other)(affString, (1, 1), affCoefs)
152 # Compute the quadratic part of the product.
153 quadPart = {(v, w): self._linear_coefs[v].T * other._linear_coefs[w]
154 for v in self._linear_coefs for w in other._linear_coefs}
156 # Don't create quadratic expressions without a quadratic part.
157 if not any(quadPart.values()):
158 affPart._symbStr = string
159 return affPart
161 # Remember a factorization into two real scalars if applicable.
162 # NOTE: If the user enters a multiplication a*b of two scalar affine
163 # expressions, then we have, at this point, self == a.T == a
164 # and other == b.conj.conj == b.
165 if len(self) == 1 and len(other) == 1 \
166 and self.isreal and other.isreal:
167 factors = (self.refined, other.refined)
168 else:
169 factors = None
171 return QuadraticExpression(
172 string, quadPart, affPart, scalarFactors=factors)
173 else:
174 return BiaffineExpression.__or__(self, other)
176 @convert_operands(rMatMul=True)
177 @refine_operands(stop_at_affine=True)
178 def __mul__(self, other):
179 if isinstance(other, ComplexAffineExpression) \
180 and not self.constant and not other.constant:
181 # If the result is scalar, allow for quadratic terms.
182 if self._shape[0] == 1 and other._shape[1] == 1 \
183 and self._shape[1] == other._shape[0]:
184 result = self.T.__or__(other.conj)
186 # NOTE: __or__ always creates a fresh expression.
187 result._symbStr = glyphs.clever_mul(self.string, other.string)
189 return result
190 else:
191 raise NotImplementedError(
192 "PICOS does not support multidimensional quadratic "
193 "expressions at this point. More precisely, one factor must"
194 " be constant or the result must be scalar.")
195 else:
196 return BiaffineExpression.__mul__(self, other)
198 @convert_operands(sameShape=True)
199 @refine_operands(stop_at_affine=True)
200 def __xor__(self, other):
201 if isinstance(other, ComplexAffineExpression) \
202 and not self.constant and not other.constant:
203 # If the result is scalar, allow for quadratic terms.
204 if self._shape == (1, 1):
205 result = self.__or__(other.conj)
207 # NOTE: __or__ always creates a fresh expression.
208 result._symbStr = glyphs.hadamard(self.string, other.string)
210 return result
211 else:
212 raise NotImplementedError(
213 "PICOS does not support multidimensional quadratic "
214 "expressions at this point. More precisely, one factor must"
215 " be constant or the result must be scalar.")
216 else:
217 return BiaffineExpression.__xor__(self, other)
219 # TODO: Create a quadratic expression from a scalar Kronecker prod.
221 # --------------------------------------------------------------------------
222 # Method overridings for BiaffineExpression: Unary operators.
223 # --------------------------------------------------------------------------
225 @cached_property
226 def real(self):
227 """Override :meth:`~.exp_biaffine.BiaffineExpression.real`.
229 The result is returned as an :meth:`AffineExpression`.
230 """
231 return AffineExpression(glyphs.real(self.string), self._shape,
232 {vars: coef.real() for vars, coef in self._coefs.items()})
234 @cached_property
235 def imag(self):
236 """Override :meth:`~.exp_biaffine.BiaffineExpression.imag`.
238 The result is returned as an :meth:`AffineExpression`.
239 """
240 return AffineExpression(glyphs.imag(self.string), self._shape,
241 {vars: coef.imag() for vars, coef in self._coefs.items()})
243 # --------------------------------------------------------------------------
244 # Additional unary operators.
245 # --------------------------------------------------------------------------
247 @cached_unary_operator
248 def __abs__(self):
249 from . import Norm
251 return Norm(self)
253 # --------------------------------------------------------------------------
254 # Constraint-creating operators, and _predict.
255 # --------------------------------------------------------------------------
257 @classmethod
258 def _predict(cls, subtype, relation, other):
259 assert isinstance(subtype, cls.Subtype)
261 from .set import Set
263 if relation == operator.__eq__:
264 if issubclass(other.clstype, ComplexAffineExpression):
265 return ComplexAffineConstraint.make_type(dim=subtype.dim)
266 elif relation == operator.__lshift__:
267 if issubclass(other.clstype, ComplexAffineExpression):
268 return ComplexLMIConstraint.make_type(int(subtype.dim**0.5))
269 elif issubclass(other.clstype, Set):
270 other_type = ExpressionType(cls, subtype)
271 return other.predict(operator.__rshift__, other_type)
272 elif relation == operator.__rshift__:
273 if issubclass(other.clstype, ComplexAffineExpression):
274 return ComplexLMIConstraint.make_type(int(subtype.dim**0.5))
276 return NotImplemented
278 @convert_operands(sameShape=True)
279 @validate_prediction
280 @refine_operands()
281 def __eq__(self, other):
282 if isinstance(other, ComplexAffineExpression):
283 return ComplexAffineConstraint(self, other)
284 else:
285 return NotImplemented
287 # Since we define __eq__, __hash__ is not inherited. Do this manually.
288 __hash__ = Expression.__hash__
290 def _lshift_implementation(self, other):
291 if isinstance(other, ComplexAffineExpression):
292 return ComplexLMIConstraint(self, Constraint.LE, other)
293 else:
294 return NotImplemented
296 def _rshift_implementation(self, other):
297 if isinstance(other, ComplexAffineExpression):
298 return ComplexLMIConstraint(self, Constraint.GE, other)
299 else:
300 return NotImplemented
302 # --------------------------------------------------------------------------
303 # Interface for PICOS-internal use.
304 # --------------------------------------------------------------------------
306 def sparse_rows(self, varOffsetMap):
307 r"""Yield a sparse list representation of the expression.
309 This is similar to :meth:`sparse_matrix_form` (with default arguments)
310 but instead of returning :math:`A` and :math:`b` at once, this yields
311 for every row of :math:`[A \mid b]`, each representing a scalar entry of
312 the expression's vectorization, a triplet containing a list of column
313 indices and values of that row of :math:`A` and the entry of :math:`b`.
315 :param varOffsetMap:
316 Maps variables to column offsets.
318 :yields tuple(list, list, float):
319 Triples ``(J, V, c)`` where ``J`` contains column indices
320 (representing scalar variables), ``V`` contains coefficients for
321 each column index, and where ``c`` is a constant term.
322 """
323 A, b = self.sparse_matrix_form(varOffsetMap, dense_b=True)
324 R, C, V = A.T.CCS
326 for r in range(len(self)):
327 u, v = R[r:r + 2]
328 yield list(C[u:v]), list(V[u:v]), b[r]
330 def sparse_matrix_form(
331 self, varOffsetMap, *, offset=0, padding=0, dense_b=False):
332 """Return a representation suited for embedding in constraint matrices.
334 This computes a sparse matrix :math:`A` and a sparse column vector
335 :math:`b` such that :math:`Ax + b` represents the vectorized expression,
336 where :math:`x` is a vertical concatenation of a number of variables,
337 including those that appear in the expression. The size and ordering of
338 :math:`x` is given through ``varOffsetMap``, which maps PICOS variables
339 to their starting position within :math:`x`.
341 If the optional parameters ``offset`` and ``padding`` are given, then
342 both :math:`A` and :math:`b` are padded with zero rows from above and
343 below, respectively.
345 This method is used by PICOS internally to assemble constraint matrices.
347 :param dict varOffsetMap:
348 Maps variables to column offsets.
349 :param int offset:
350 Number of zero rows to insert at the top of :math:`A` and :math:`b`.
351 :param int offset:
352 Number of zero rows to insert at the bottom of :math:`A` and
353 :math:`b`.
354 :param bool dense_b:
355 Whether to return :math:`b` as a dense vector. Not compatible with
356 nonzero ``offset`` or ``padding``.
358 :returns tuple(cvxopt.spmatrix):
359 A pair ``(A, b)`` of CVXOPT sparse matrices representing the matrix
360 :math:`A` and the column vector :math:`b`. (If ``dense_b=True``,
361 then ``b`` is returned as a dense CVXOPT column vector instead.)
362 """
363 lin = self._sparse_linear_coefs
364 cst = self._constant_coef
366 tc = self._typecode
368 k = len(self)
369 m = offset + k + padding
370 n = sum(var.dim for var in varOffsetMap)
372 ordered_vars = sorted(varOffsetMap, key=varOffsetMap.__getitem__)
374 blocks = []
375 for var in ordered_vars:
376 if var in lin:
377 coef = lin[var]
378 else:
379 coef = cvxopt.spmatrix([], [], [], size=(k, var.dim), tc=tc)
381 blocks.append([coef])
383 if blocks:
384 A = cvxopt.sparse(blocks, tc=tc)
385 else:
386 A = cvxopt.spmatrix([], [], [], size=(k, 0), tc=tc)
388 b = cvxopt.matrix(cst, tc=tc) if dense_b else cvxopt.sparse(cst, tc=tc)
390 if offset or padding:
391 if dense_b:
392 raise ValueError("Refusing to return a dense vector if a "
393 "nonzero offset or padding is given.")
395 A = cvxopt.sparse(
396 [
397 cvxopt.spmatrix([], [], [], size=(offset, n), tc=tc),
398 A,
399 cvxopt.spmatrix([], [], [], size=(padding, n), tc=tc)
400 ], tc=tc
401 )
403 b = cvxopt.sparse(
404 [
405 cvxopt.spmatrix([], [], [], size=(offset, 1), tc=tc),
406 b,
407 cvxopt.spmatrix([], [], [], size=(padding, 1), tc=tc)
408 ], tc=tc
409 )
411 assert A.size == (m, n)
412 assert b.size == (m, 1)
414 return A, b
416 def scipy_sparse_matrix_form(
417 self, varOffsetMap, *, offset=0, padding=0, dense_b=False):
418 """Like :meth:`sparse_matrix_form` but returns SciPy types.
420 See :meth:`sparse_matrix_form` for details and arguments.
422 :returns tuple(scipy.sparse.csc_matrix):
423 A pair ``(A, b)`` of SciPy sparse matrices in CSC format
424 representing the matrix :math:`A` and the column vector :math:`b`.
425 (If ``dense_b=True``, then ``b`` is returned as a 1-D :class:`NumPy
426 array <numpy.ndarray>` instead.)
428 :raises ModuleNotFoundError:
429 If the optional dependency :mod:`scipy` is not installed.
430 """
431 import scipy.sparse
433 lin = self._linear_coefs
434 cst = self._constant_coef
436 dtype = type(cst[0])
438 k = len(self)
439 m = offset + k + padding
440 n = sum(var.dim for var in varOffsetMap)
442 ordered_vars = sorted(varOffsetMap, key=varOffsetMap.__getitem__)
444 blocks = []
445 for var in ordered_vars:
446 if var in lin:
447 coef = cvx2csc(lin[var])
448 else:
449 coef = scipy.sparse.csc_matrix((k, var.dim), dtype=dtype)
451 blocks.append(coef)
453 if blocks:
454 A = scipy.sparse.hstack(blocks, format="csc", dtype=dtype)
455 else:
456 A = scipy.sparse.csc_matrix((k, 0), dtype=dtype)
458 b = numpy.ravel(cvx2np(cst)) if dense_b else cvx2csc(cst)
460 if offset or padding:
461 if dense_b:
462 raise ValueError("Refusing to return a dense vector if a "
463 "nonzero offset or padding is given.")
465 A = scipy.sparse.vstack(
466 [
467 scipy.sparse.csc_matrix((offset, n), dtype=dtype),
468 A,
469 scipy.sparse.csc_matrix((padding, n), dtype=dtype)
470 ], format="csc", dtype=dtype
471 )
473 b = scipy.sparse.vstack(
474 [
475 scipy.sparse.csc_matrix((offset, 1), dtype=dtype),
476 b,
477 scipy.sparse.csc_matrix((padding, 1), dtype=dtype)
478 ], format="csc", dtype=dtype
479 )
481 assert A.shape == (m, n)
482 assert b.shape in ((m,), (m, 1))
484 return A, b
487class AffineExpression(ComplexAffineExpression):
488 """A multidimensional real affine expression."""
490 # --------------------------------------------------------------------------
491 # Method overridings for BiaffineExpression.
492 # --------------------------------------------------------------------------
494 @property
495 def isreal(self):
496 """Always true for :class:`AffineExpression` instances.""" # noqa
497 return True
499 @property
500 def real(self):
501 """The :class:`AffineExpression` as is.""" # noqa
502 return self
504 @cached_property
505 def imag(self):
506 """A zero of same shape as the :class:`AffineExpression`.""" # noqa
507 return self._basetype.zero(self._shape)
509 @property
510 def conj(self):
511 """The :class:`AffineExpression` as is.""" # noqa
512 return self
514 @property
515 def H(self):
516 """The regular transpose of the :class:`AffineExpression`.""" # noqa
517 return self.T
519 # --------------------------------------------------------------------------
520 # Method overridings for ComplexAffineExpression.
521 # --------------------------------------------------------------------------
523 @classmethod
524 def _get_basetype(cls):
525 return AffineExpression
527 @classmethod
528 def _get_typecode(cls):
529 return "d"
531 def _get_refined(self):
532 return self
534 @convert_operands(sameShape=True, allowNone=True)
535 def _set_value(self, value):
536 if value is None:
537 for var in self._linear_coefs:
538 var.value = None
539 return
541 if not isinstance(value, AffineExpression) or not value.constant:
542 raise TypeError("Cannot set the value of {} to {}: Not real or not "
543 "a constant.".format(repr(self), repr(value)))
545 if self.constant:
546 raise TypeError("Cannot set the value on a constant expression.")
548 y = cvx2np(value._constant_coef)
550 A = []
551 for var, coef in self._linear_coefs.items():
552 A.append(cvx2np(coef))
553 assert A
555 A = numpy.hstack(A)
556 b = y - cvx2np(self._constant_coef)
558 try:
559 solution, residual, _, _ = numpy.linalg.lstsq(A, b, rcond=None)
560 except numpy.linalg.LinAlgError as error:
561 raise RuntimeError("Setting a value on {} by means of a least-"
562 "squares solution failed.".format(self.string)) from error
564 if not numpy.allclose(residual, 0):
565 raise ValueError("Setting a value on {} failed: No exact solution "
566 "to the associated linear system found.".format(self.string))
568 offset = 0
569 for var in self._linear_coefs:
570 var.internal_value = solution[offset:offset+var.dim]
571 offset += var.dim
573 # --------------------------------------------------------------------------
574 # Additional unary operators.
575 # --------------------------------------------------------------------------
577 @cached_property
578 def exp(self):
579 """The exponential function applied to the expression.""" # noqa
580 from . import SumExponentials
581 return SumExponentials(self)
583 @cached_property
584 def log(self):
585 """The Logarithm of the expression.""" # noqa
586 from . import Logarithm
587 return Logarithm(self)
589 # --------------------------------------------------------------------------
590 # Constraint-creating operators, and _predict.
591 # --------------------------------------------------------------------------
593 @classmethod
594 def _predict(cls, subtype, relation, other):
595 assert isinstance(subtype, cls.Subtype)
597 if relation in (operator.__eq__, operator.__le__, operator.__ge__):
598 if issubclass(other.clstype, AffineExpression):
599 return AffineConstraint.make_type(
600 dim=subtype.dim, eq=(relation is operator.__eq__))
601 elif relation == operator.__lshift__:
602 if issubclass(other.clstype, AffineExpression):
603 return LMIConstraint.make_type(int(subtype.dim**0.5))
604 elif issubclass(other.clstype, ComplexAffineExpression):
605 return ComplexLMIConstraint.make_type(int(subtype.dim**0.5))
606 elif relation == operator.__rshift__:
607 if issubclass(other.clstype, AffineExpression):
608 return LMIConstraint.make_type(int(subtype.dim**0.5))
609 elif issubclass(other.clstype, ComplexAffineExpression):
610 return ComplexLMIConstraint.make_type(int(subtype.dim**0.5))
612 return NotImplemented
614 @convert_operands(sameShape=True)
615 @validate_prediction
616 @refine_operands()
617 def __le__(self, other):
618 if isinstance(other, AffineExpression):
619 return AffineConstraint(self, Constraint.LE, other)
620 else:
621 return NotImplemented
623 @convert_operands(sameShape=True)
624 @validate_prediction
625 @refine_operands()
626 def __ge__(self, other):
627 if isinstance(other, AffineExpression):
628 return AffineConstraint(self, Constraint.GE, other)
629 else:
630 return NotImplemented
632 @convert_operands(sameShape=True)
633 @validate_prediction
634 @refine_operands()
635 def __eq__(self, other):
636 if isinstance(other, AffineExpression):
637 return AffineConstraint(self, Constraint.EQ, other)
638 else:
639 return NotImplemented
641 # Since we define __eq__, __hash__ is not inherited. Do this manually.
642 __hash__ = Expression.__hash__
644 def _lshift_implementation(self, other):
645 if isinstance(other, AffineExpression):
646 return LMIConstraint(self, Constraint.LE, other)
647 elif isinstance(other, ComplexAffineExpression):
648 return ComplexLMIConstraint(self, Constraint.LE, other)
649 else:
650 return NotImplemented
652 def _rshift_implementation(self, other):
653 if isinstance(other, AffineExpression):
654 return LMIConstraint(self, Constraint.GE, other)
655 elif isinstance(other, ComplexAffineExpression):
656 return ComplexLMIConstraint(self, Constraint.GE, other)
657 else:
658 return NotImplemented
661def Constant(name_or_value, value=None, shape=None):
662 """Create a constant PICOS expression.
664 Loads the given numeric value as a constant
665 :class:`~picos.expressions.ComplexAffineExpression` or
666 :class:`~picos.expressions.AffineExpression`, depending on the value.
667 Optionally, the value is broadcasted or reshaped according to the shape
668 argument.
670 :param str name_or_value: Symbolic string description of the constant. If
671 :obj:`None` or the empty string, a string will be generated. If this is
672 the only positional parameter (i.e.``value`` is not given), then this
673 position is used as the value argument instead!
674 :param value: The numeric constant to load.
676 See :func:`~.data.load_data` for supported data formats and broadcasting and
677 reshaping rules.
679 :Example:
681 >>> from picos import Constant
682 >>> Constant(1)
683 <1×1 Real Constant: 1>
684 >>> Constant(1, shape=(2, 2))
685 <2×2 Real Constant: [1]>
686 >>> Constant("one", 1)
687 <1×1 Real Constant: one>
688 >>> Constant("J", 1, (2, 2))
689 <2×2 Real Constant: J>
690 """
691 if value is None:
692 value = name_or_value
693 name = None
694 else:
695 name = name_or_value
697 value, valStr = load_data(value, shape)
699 if value.typecode == "z":
700 cls = ComplexAffineExpression
701 else:
702 cls = AffineExpression
704 return cls(name if name else valStr, value.size, {(): value})
707# --------------------------------------
708__all__ = api_end(_API_START, globals())