Coverage for picos/expressions/data.py: 70.50%
827 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"""Functions to load and work with raw numeric data."""
22import functools
23import inspect
24import math
25import re
26import sys
27from fractions import Fraction
28from functools import lru_cache
30import cvxopt
31import numpy
33from .. import glyphs, settings
34from ..apidoc import api_end, api_start
36_API_START = api_start(globals())
37# -------------------------------
40#: Maximum entrywise absolute deviation allowed for numeric equality checks.
41TOLERANCE = 1e-6
44def load_shape(shape, squareMatrix=False, wildcards=False):
45 """Parse the argument as a matrix shape.
47 PICOS uses this function whenever you supply a shape parameter to a method.
49 A scalar argument is treated as the length of a column-vector. If the shape
50 contains :obj:`None`, it is treated as a wildcard (any dimensionality).
52 :param bool squareMatrix: If :obj:`True`, a scalar argument is treated as
53 the side/diagonal length of a square matrix, and any other argument is
54 validated to be square. If :obj:`False`, a scalar argument is treated
55 as the length of a column vector.
56 :param bool wildcards: Whether the wildcard token :obj:`None` is allowed.
57 """
58 if shape is None:
59 shape = (None, None)
60 elif isinstance(shape, int):
61 if squareMatrix:
62 shape = (shape, shape)
63 else:
64 shape = (shape, 1)
65 elif not isinstance(shape, tuple) and not isinstance(shape, list):
66 raise TypeError("Shapes must be given as None, int, tuple or list.")
67 elif len(shape) == 1:
68 shape = (shape[0], 1)
69 elif len(shape) == 0:
70 shape = (1, 1)
71 elif len(shape) != 2:
72 raise TypeError("Shapes must be two-dimensional.")
74 shape = (
75 None if shape[0] is None else int(shape[0]),
76 None if shape[1] is None else int(shape[1]))
78 if not wildcards and None in shape:
79 raise ValueError("Invalid shape (wildcards not allowed): {}."
80 .format(glyphs.shape(shape)))
82 if 0 in shape:
83 raise ValueError("Invalid shape (zero-dimensional axis): {}."
84 .format(glyphs.shape(shape)))
86 if squareMatrix and shape[0] != shape[1]:
87 raise ValueError("Invalid shape for a square matrix: {}."
88 .format(glyphs.shape(shape)))
90 return shape
93def blend_shapes(baseShape, defaultShape):
94 """Replace wildcards in one shape with entries of the other.
96 :param baseShape: Primary shape, usually with wildcards.
97 :type baseShape: tuple(int or None)
98 :param defaultShape: Secondary shape with fallback entries.
99 :type defaultShape: tuple(int or None)
100 """
101 return (
102 defaultShape[0] if baseShape[0] is None else baseShape[0],
103 defaultShape[1] if baseShape[1] is None else baseShape[1])
106def should_be_sparse(shape, numNonZero):
107 """Decide whether a matrix is considered sparse.
109 :param tuple(int) shape: The shape of the matrix in question.
110 :param int numNonZero: Number of non-zero elements of the matrix.
111 """
112 n, m = shape
113 l = (n*m)**0.5
114 return numNonZero < l * math.log(l)
117_LOAD_DATA_REGEX = re.compile("^{}{}{}$".format(
118 r"([\-0-9.+j]+)?", # Leading coefficient.
119 "(" + "|".join(( # Matrix type:
120 r"e_[0-9]+(?:,[0-9]+)?", # Single nonzero element.
121 r"\|[\-0-9.+j]+\|", # All equal elements.
122 "I", # Identity matrix.
123 )) + ")",
124 r"(\([0-9]+(?:,[0-9]+)?\))?" # Fixed shape.
125))
128def is_scipy_spmat(value):
129 """Report whether value is a SciPy sparse matrix without importing scipy."""
130 return hasattr(value, "__module__") \
131 and value.__module__.startswith("scipy.sparse")
134def load_data(value, shape=None, typecode=None, sparse=None,
135 alwaysCopy=True, legacy=False):
136 r"""Load a constant numeric data value as a CVXOPT (sparse) matrix.
138 As a user, you never need to call this manually, but you should be aware
139 that PICOS uses this function on any raw data you supply as an operand when
140 working with PICOS expressions. For instance, you can just add a NumPy
141 matrix or an algebraic string such as ``"I"`` to such an expression without
142 worrying about any conversion.
144 :Supported data types:
146 - A NumPy matrix: :class:`numpy.ndarray` or :class:`numpy.matrix`.
147 - A SciPy sparse matrix: All from :mod:`scipy.sparse`.
148 - A CVXOPT matrix: :obj:`cvxopt.matrix` or :obj:`cvxopt.spmatrix`.
149 - A constant PICOS expression: :class:`~.exp_affine.AffineExpression` or
150 :class:`~.exp_affine.ComplexAffineExpression`.
151 - A Python scalar: :class:`int`, :class:`float` or :class:`complex`.
152 - A flat :class:`tuple` or :class:`list` containing scalars or a
153 :class:`range`, all representing a column vector.
154 - A nested :class:`tuple` or :class:`list` containing scalars. The outer
155 level represents rows and the inner level represents the rows'
156 entries. Allows you to define a :math:`2 \times 3` matrix like this:
158 .. code-block:: python
160 A = [[1, 2, 3],
161 [4, 5, 6]]
163 - A verbatim string description, with rows separated by newline and
164 columns separated by whitespace. The same :math:`2 \times 3` matrix as
165 above:
167 .. code-block:: python
169 A = '''1 2 3
170 4 5 6'''
172 - An algebraic string description:
174 .. list-table::
175 :widths: 1 99
177 * - ``"|a|"``
178 - A matrix with all entries equal to :math:`a`.
179 * - ``"|a|(m,n)"``
180 - A :math:`m \times n` matrix with all entries equal to :math:`a`.
181 * - ``"e_i,j(m,n)"``
182 - A :math:`m \times n` matrix with a :math:`1` at :math:`(i,j)`,
183 indexed from :math:`(1,1)`` to :math:`(m,n)``.
184 * - ``"e_i(m,n)"``
185 - A :math:`m \times n` matrix with a single :math:`1` on the
186 :math:`i`-th coordinate, indexed from :math:`1` in column-major
187 order.
188 * - ``"I"``
189 - The identity matrix.
190 * - ``"I(n)"``
191 - The :math:`n \times n` identiy matrix.
192 * - ``"a…"``
193 - The matrix given by ``"…"`` but multiplied by :math:`a`.
195 Different matrix operations such as addition or multiplication have
196 different requirements with respect to the operands' shapes. The shape
197 of any PICOS expression involved will always be maintained. But when an
198 operation involves both a PICOS expression and raw data, then PICOS will
199 try to broadcast or reshape the raw data such that the operation can be
200 performed.
202 :Broadcasting and reshaping rules:
204 - An input vector without a second axis (for instance a non-nested
205 :class:`tuple` or :class:`list` or a :class:`range`) is interpreted as
206 a row vector if the target shape is of the form ``(None, n)`` with
207 :math:`n > 1`, otherwise it is interpreted as a column vector.
208 - If the target shape is :obj:`None` or ``(None, None)``, then the
209 input's shape is maintained.
210 - If the target shape contains :obj:`None` exactly once, that occurence
211 is replaced by the respective dimensionality of the input data shape.
212 - A scalar is copied (broadcasted) to fill the target shape.
213 - A column (row) vector is copied horizontally (vertically) if its
214 length matches the target row (column) count.
215 - Reshaping from matrix to vector: A matrix is vectorized in
216 column-major (row-major) order if the target shape is a column (row)
217 vector whose length matches the number of matrix entries.
218 - Reshaping from vector to matrix: A column (row) vector is treated as
219 the column-major (row-major) vectorization of the target matrix if the
220 former's length matches the number of the latter's entries.
221 - All other combinations of input and target shape raise an exception.
222 - When an algebraic string description specifies no shape, then the
223 shape argument must be supplied. When both the string and the shape
224 argument specify a shape, then they must be consistent (no
225 broadcasting or reshaping is applied in this case).
227 :param shape: The shape of the resulting matrix. If the input data is of
228 another shape, broadcasting or reshaping is used if possible, otherwise
229 an exception is raised. An integer is treated as the length of a column
230 vector. If this is :obj:`None`, then the target's shape is the input's.
231 If only the target number of rows or columns is :obj:`None`, then only
232 that quantity is chosen according to the input.
233 :type shape: int or tuple or list or None
235 :param str typecode: The numeric type of the resulting matrix. Either
236 ``'d'`` (float), ``'z'`` (complex) or ``'i'`` (integer). If the input
237 data is not already of this type, then it will be converted if possible.
238 If this is not possible, then an exception is raised. If this is
239 :obj:`None`, then the output type is chosen based on the input data.
241 :param sparse: If :obj:`True`, a sparse matrix is returned. If :obj:`False`,
242 a dense matrix is returned. If :obj:`None`, it depends on the sparsity
243 pattern of the input data. If the typecode argument is ``'i'``, then
244 a value of :obj:`True` is not allowed and the returned matrix is dense.
245 :type sparse: bool or None
247 :param bool alwaysCopy: If :obj:`True`, then a copy of the input data is
248 returned even if it already equals the output data. If :obj:`False`,
249 the input value can be returned as such if it is already a CVXOPT matrix
250 of the desired shape, typecode, and sparsity.
252 :param bool legacy: Be compatible with the old ``retrieve_matrix`` function.
253 In particular, if the target shape contains :obj:`None` exactly once and
254 the input data is scalar, treat this as a matrix multiplication case and
255 return the scalar times an identity matrix of appropriate size.
257 :returns: A :class:`tuple` whose first entry is the loaded matrix and whose
258 second argument is a short string for representing the data within
259 algebraic expression strings.
261 :Example:
263 >>> from picos.expressions.data import load_data
264 >>> # Data as (nested) list:
265 >>> load_data([1,2,3])
266 (<3x1 matrix, tc='i'>, '[3×1]')
267 >>> load_data([[1,2,3]])
268 (<1x3 matrix, tc='i'>, '[1×3]')
269 >>> A = [[1,2,3],
270 ... [4,5,6]]
271 >>> load_data(A)
272 (<2x3 matrix, tc='i'>, '[2×3]')
273 >>> # Data as string:
274 >>> value, string = load_data('e_14(7,2)')
275 >>> print(string)
276 [7×2:e_7,2]
277 >>> print(value) #doctest: +NORMALIZE_WHITESPACE
278 [ 0 0 ]
279 [ 0 0 ]
280 [ 0 0 ]
281 [ 0 0 ]
282 [ 0 0 ]
283 [ 0 0 ]
284 [ 0 1.00e+00]
285 >>> load_data('5.3I', (2,2))
286 (<2x2 sparse matrix, tc='d', nnz=2>, '5.3·I')
287 """
288 from .exp_affine import ComplexAffineExpression
289 from .expression import Expression
291 def load_sparse(V, I, J, shape, typecode):
292 """Create a CVXOPT sparse matrix."""
293 # HACK: Work around CVXOPT not supporting integer sparse matrices:
294 # Create a real sparse matrix for now (will be converted later).
295 # Note that users may not request both sparsity and integrality.
296 typecode = "d" if typecode == "i" else typecode
298 try:
299 if typecode:
300 return cvxopt.spmatrix(V, I, J, shape, typecode)
301 else:
302 return cvxopt.spmatrix(V, I, J, shape)
303 except TypeError as error:
304 # Attempt to convert complex typed but real valued input to a real
305 # typed output matrix.
306 if typecode == "d":
307 realV = [x.real for x in V]
308 if realV == V:
309 try:
310 return cvxopt.spmatrix(realV, I, J, shape, typecode)
311 except Exception:
312 pass
314 raise TypeError(
315 "Failed to create a CVXOPT sparse matrix of shape {} and type "
316 "'{}'.".format(glyphs.shape(shape), typecode)) from error
318 def load_dense(value, shape, typecode):
319 """Create a CVXOPT dense matrix."""
320 try:
321 if typecode:
322 return cvxopt.matrix(value, shape, typecode)
323 else:
324 return cvxopt.matrix(value, shape)
325 except TypeError as error:
326 # Attempt to convert complex (real) typed but real/integer (integer)
327 # valued input to a real/integer (integer) typed output matrix.
328 if typecode in "id":
329 try:
330 complexValue = cvxopt.matrix(value, shape, "z")
331 except Exception:
332 pass
333 else:
334 if not any(complexValue.imag()):
335 if typecode == "d":
336 return complexValue.real()
337 else:
338 realData = list(complexValue.real())
339 intData = [int(x) for x in realData]
340 if intData == realData:
341 try:
342 return cvxopt.matrix(intData, shape, "i")
343 except Exception:
344 pass
346 raise TypeError(
347 "Failed to create a CVXOPT dense matrix of shape {} and type "
348 "'{}'.".format(glyphs.shape(shape), typecode)) from error
350 def simple_vector_as_row(shape):
351 """Whether a single-axis vector should be a row or column vector."""
352 return shape[0] is None and shape[1] is not None and shape[1] > 1
354 def broadcast_error():
355 """Raise a broadcasting related :class:TypeError."""
356 raise TypeError("Cannot broadcast or reshape from {} to {}{}."
357 .format(glyphs.shape(inShape), glyphs.shape(shape), " read as {}"
358 .format(glyphs.shape(outShape)) if shape != outShape else ""))
360 def scalar(x, typecode=None):
361 """Load a scalar as either a complex, float, or int."""
362 x = complex(x)
364 if typecode == "z":
365 return x # We are content with complex.
366 else:
367 if not x.imag:
368 x = x.real
370 if typecode == "d":
371 return x # We are content with real.
372 else:
373 if int(x) == x:
374 return int(x)
375 elif not typecode:
376 return x # We are not strict and real is the best.
377 elif not typecode:
378 return x # We are not strict and complex is the best.
380 raise TypeError(
381 "Cannot cast {} according to typecode '{}'.".format(x, typecode))
383 # Normalize the shape argument to a two-dimensional tuple of int or None,
384 # where None means any dimensionality.
385 shape = load_shape(shape, wildcards=True)
387 # Validate the typecode argument.
388 if typecode is not None and typecode not in "dzi":
389 raise ValueError("Typecode argument not a valid CVXOPT typecode: {}."
390 .format(typecode))
392 # Validate the sparsity argument.
393 if sparse not in (None, True, False):
394 raise ValueError("Sparsity argument must be True, False or None.")
396 # CVXOPT sparse matrices may not be integer typed.
397 if typecode == "i":
398 if sparse:
399 raise TypeError("Sparse integer matrices are not implemented with "
400 "CVXOPT, which PICOS uses as its matrix backed.")
401 else:
402 sparse = False
404 # Conversions must retrieve the input shape for further processing.
405 inShape = None
407 # Allow conversions to specify their own string description.
408 string = None
410 # Convert from range to list.
411 if isinstance(value, range):
412 value = list(value)
414 # Try to refine a PICOS expression to a constant affine expression.
415 if isinstance(value, Expression):
416 value = value.refined
418 # Convert from a SciPy sparse to a CVXOPT sparse matrix.
419 if is_scipy_spmat(value):
420 value = sp2cvx(value)
421 alwaysCopy = False
423 # Convert the data to a CVXOPT (sparse) matrix of proper shape and type.
424 if isinstance(value, (int, float, complex, numpy.number)):
425 # Convert NumPy numbers to Python ones.
426 if isinstance(value, numpy.number):
427 if isinstance(value, numpy.integer):
428 value = int(value)
429 elif isinstance(value, numpy.floating):
430 value = float(value)
431 elif isinstance(value, numpy.complexfloating):
432 value = complex(value)
433 else:
434 assert False, "Unexpected NumPy numeric type {}.".format(
435 type(value).__name__)
437 # CVXOPT is limited by the system's word length.
438 if isinstance(value, int):
439 if value > sys.maxsize or value < -sys.maxsize - 1:
440 raise ValueError("The number {} is too large to be loaded by "
441 "PICOS.".format(value))
443 inShape = (1, 1)
444 outShape = blend_shapes(shape, inShape)
446 string = glyphs.scalar(value)
447 if value and outShape != (1, 1): # Don't use the matrix glyph on zero.
448 string = glyphs.matrix(string)
450 if not value and sparse is not False:
451 value = load_sparse([], [], [], outShape, typecode)
452 else:
453 value = load_dense(value, outShape, typecode)
455 if sparse:
456 value = cvxopt.sparse(value)
457 elif isinstance(value, cvxopt.matrix):
458 # NOTE: Since the input is already a CVXOPT data type, it is possible
459 # that it can be returned without modification. The 'alwaysCopy'
460 # parameter, if set to True, prevents this. As soon as a copy is
461 # made during processing of the input we set 'alwaysCopy' to False
462 # to prevent an unnecessary second copy of the processed input to
463 # be made later on.
465 # If sparse is requested, let CVXOPT handle the transformation first.
466 if sparse:
467 value = cvxopt.sparse(value)
468 return load_data(value, shape, typecode, sparse, False, legacy)
470 # Refine the output shape.
471 inShape = value.size
472 outShape = blend_shapes(shape, inShape)
474 # Define shorthands.
475 inLength = inShape[0] * inShape[1]
476 outLength = outShape[0] * outShape[1]
478 # If the input is stored as complex and no complex output is explicitly
479 # requested, try to remove an imaginary part of zero.
480 if value.typecode == "z" and typecode != "z" and not any(value.imag()):
481 value = value.real()
482 alwaysCopy = False
484 # Broadcast or reshape the data if necessary and possible.
485 reshaped = True
486 if inShape == outShape:
487 reshaped = False
488 elif inShape == (1, 1):
489 # Broadcast a scalar.
490 value = load_dense(value[0], outShape, typecode)
491 elif inShape == (outShape[0], 1):
492 # Copy columns horizontally.
493 value = load_dense(list(value)*outShape[1], outShape, typecode)
494 elif inShape == (1, outShape[1]):
495 # Copy rows vertically.
496 value = load_dense(
497 list(value)*outShape[0], (outShape[1], outShape[0]), typecode).T
498 elif (inLength, 1) == outShape:
499 # Vectorize in column-major order.
500 if not typecode or value.typecode == typecode:
501 # Quick method.
502 value = value[:]
503 else:
504 # With typecode change.
505 value = load_dense(value, outShape, typecode)
506 elif (1, inLength) == outShape:
507 # Vectorize in row-major order.
508 if not typecode or value.typecode == typecode:
509 # Quick method.
510 value = value.T[:].T
511 else:
512 # With typecode change.
513 value = load_dense(value.T, outShape, typecode)
514 elif (outLength, 1) == inShape:
515 # Devectorize in column-major order.
516 value = load_dense(value, outShape, typecode)
517 elif (1, outLength) == inShape:
518 # Devectorize in row-major order.
519 value = load_dense(value, (outShape[1], outShape[0]), typecode).T
520 else:
521 broadcast_error()
522 if reshaped:
523 alwaysCopy = False
525 # The data now has the desired shape.
526 assert value.size == outShape
528 # Ensure the proper typecode.
529 if typecode and value.typecode != typecode:
530 value = load_dense(value, outShape, typecode)
531 alwaysCopy = False
533 # Copy the data if requested and not already a copy.
534 if alwaysCopy:
535 value = load_dense(value, outShape, typecode)
536 elif isinstance(value, cvxopt.spmatrix):
537 # NOTE: See case above.
539 # Refine the output shape.
540 inShape = value.size
541 outShape = blend_shapes(shape, inShape)
543 # Define shorthands.
544 inLength = inShape[0] * inShape[1]
545 outLength = outShape[0] * outShape[1]
547 # If the input is stored as complex and no complex output is explicitly
548 # requested, try to remove an imaginary part of zero.
549 if value.typecode == "z" and typecode != "z" and not any(value.imag()):
550 value = value.real()
551 alwaysCopy = False
553 # Broadcast or reshape the data if necessary and possible.
554 reshaped = True
555 if inShape == outShape:
556 reshaped = False
557 elif inShape == (1, 1):
558 # Broadcast a scalar.
559 if value[0] != 0:
560 value = load_dense(value[0], outShape, typecode)
561 else:
562 value = load_sparse([], [], [], outShape, typecode)
563 elif inShape == (outShape[0], 1):
564 # Copy columns horizontally.
565 V = list(value.V)*outShape[1]
566 I = list(value.I)*outShape[1]
567 J = [j for j in range(outShape[1]) for _ in range(len(value))]
568 value = load_sparse(V, I, J, outShape, typecode)
569 elif inShape == (1, outShape[1]):
570 # Copy rows vertically.
571 V = list(value.V)*outShape[0]
572 I = [i for i in range(outShape[0]) for _ in range(len(value))]
573 J = list(value.J)*outShape[0]
574 value = load_sparse(V, I, J, outShape, typecode)
575 elif (inLength, 1) == outShape:
576 # Vectorize in column-major order.
577 if not typecode or value.typecode == typecode:
578 # Quick method.
579 value = value[:]
580 else:
581 # With typecode change.
582 n, nnz = inShape[0], len(value)
583 I, J = value.I, value.J
584 I = [I[k] % n + J[k]*n for k in range(nnz)]
585 J = [0]*nnz
586 value = load_sparse(value.V, I, J, outShape, typecode)
587 elif (1, inLength) == outShape:
588 # Vectorize in row-major order.
589 if not typecode or value.typecode == typecode:
590 # Quick method.
591 value = value.T[:].T
592 else:
593 # With typecode change.
594 m, nnz = inShape[1], len(value)
595 I, J = value.I, value.J
596 J = [I[k]*m + J[k] % m for k in range(nnz)]
597 I = [0]*nnz
598 value = load_sparse(value.V, I, J, outShape, typecode)
599 elif (outLength, 1) == inShape:
600 # Devectorize in column-major order.
601 # TODO: Logic for also changing the typecode.
602 value = value[:]
603 value.size = outShape
604 elif (1, outLength) == inShape:
605 # Devectorize in row-major order.
606 # TODO: Logic for also changing the typecode.
607 value = value[:]
608 value.size = (outShape[1], outShape[0])
609 value = value.T
610 else:
611 broadcast_error()
612 if reshaped:
613 alwaysCopy = False
615 # The data now has the desired shape.
616 assert value.size == outShape
618 # Ensure the proper typecode.
619 if typecode and value.typecode != typecode:
620 # NOTE: In the case of intential loading as a dense matrix, the
621 # typecode is already set properly.
622 assert isinstance(value, cvxopt.spmatrix)
623 value = load_sparse(value.V, value.I, value.J, outShape, typecode)
624 alwaysCopy = False
626 # Return either a (sparse) copy or the input data itself.
627 if isinstance(value, cvxopt.matrix):
628 # The data was intentionally copied to a dense matrix (through
629 # broadcasting), so only convert it back to sparse if requested.
630 assert reshaped
631 if sparse:
632 value = cvxopt.sparse(value)
633 elif sparse is False:
634 value = load_dense(value, outShape, typecode)
635 elif alwaysCopy:
636 value = load_sparse(value.V, value.I, value.J, outShape, typecode)
637 elif isinstance(value, numpy.ndarray):
638 # NumPy arrays can be tensors, we don't support those.
639 if len(value.shape) > 2:
640 raise TypeError("PICOS does not support tensor data.")
642 # Refine the output shape.
643 inShape = load_shape(value.shape)
644 outShape = blend_shapes(shape, inShape)
646 # Define shorthands.
647 inLength = inShape[0] * inShape[1]
648 outLength = outShape[0] * outShape[1]
650 # If the input is one-dimensional, turn it into a column or row vector.
651 if inShape != value.shape:
652 if simple_vector_as_row(shape):
653 value = value.reshape((1, inLength))
654 else:
655 value = value.reshape(inShape)
657 # If the input is stored as complex and no complex output is explicitly
658 # requested, try to remove an imaginary part of zero.
659 if value.dtype.kind == "c" and typecode != "z" \
660 and not numpy.iscomplex(value).any():
661 value = value.real
663 # Broadcast or reshape the data if necessary and possible.
664 if inShape == outShape:
665 pass
666 elif inShape == (1, 1):
667 # Broadcast a scalar.
668 value = numpy.full(outShape, value)
669 elif inShape == (outShape[0], 1):
670 # Copy columns horizontally.
671 value = value.repeat(outShape[1], 1)
672 elif inShape == (1, outShape[1]):
673 # Copy rows vertically.
674 value = value.repeat(outShape[0], 0)
675 elif (inLength, 1) == outShape:
676 # Vectorize in column-major order.
677 value = value.reshape(outShape, order="F")
678 elif (1, inLength) == outShape:
679 # Vectorize in row-major order.
680 value = value.reshape(outShape, order="C")
681 elif (outLength, 1) == inShape:
682 # Devectorize in column-major order.
683 value = value.reshape(outShape, order="F")
684 elif (1, outLength) == inShape:
685 # Devectorize in row-major order.
686 value = value.reshape(outShape, order="C")
687 else:
688 broadcast_error()
690 # The data now has the desired shape.
691 assert value.shape == outShape
693 # Decide on whether to create a dense or sparse matrix.
694 outSparse = should_be_sparse(outShape, numpy.count_nonzero(value)) \
695 if sparse is None else sparse
697 # Convert to CVXOPT.
698 if outSparse:
699 I, J = value.nonzero()
700 V = value[I, J]
701 value = load_sparse(V, I, J, outShape, typecode)
702 else:
703 value = load_dense(value, outShape, typecode)
704 elif isinstance(value, ComplexAffineExpression):
705 # Must be constant.
706 if not value.constant:
707 raise ValueError("Cannot load the nonconstant expression {} as a "
708 "constant data value.".format(value.string))
710 # Retrieve a copy of the numeric value.
711 value = value.value
713 # NOTE: alwaysCopy=False as the value is already a copy.
714 return load_data(value, shape, typecode, sparse, False, legacy)
715 elif isinstance(value, str) and re.search(r"\s", value):
716 value = value.strip().splitlines()
717 value = [row.strip().split() for row in value]
718 value = [[scalar(x, typecode) for x in row] for row in value]
720 return load_data(value, shape, typecode, sparse, alwaysCopy, legacy)
721 elif isinstance(value, str):
722 # Verbatim matrix strings with only one element fall throught to this
723 # case as they don't (have to) contain a whitespace character.
724 try:
725 value = scalar(value, typecode)
726 except ValueError:
727 pass
728 else:
729 return load_data(value, shape, typecode, sparse, alwaysCopy, legacy)
731 match = _LOAD_DATA_REGEX.match(value)
733 if not match:
734 raise ValueError("The string '{}' could not be parsed as a matrix."
735 .format(value))
737 # Retrieve the string tokens.
738 tokens = match.groups()
739 assert len(tokens) == 3
740 factor, base, inShape = tokens
741 assert base is not None
743 # Convert the factor.
744 factor = scalar(factor) if factor else 1
746 # Determine whether the matrix will be square.
747 square = base == "I"
749 # Convert the shape.
750 if inShape:
751 inShape = inShape[1:-1].split(",")
752 if len(inShape) == 1:
753 inShape *= 2
754 inShape = (int(inShape[0]), int(inShape[1]))
756 if blend_shapes(shape, inShape) != inShape:
757 raise ValueError(
758 "Inconsistent shapes for matrix given as '{}' with expected"
759 " shape {}.".format(value, glyphs.shape(shape)))
761 outShape = inShape
762 elif None in shape:
763 if square and shape != (None, None):
764 outShape = blend_shapes(shape, (shape[1], shape[0]))
765 assert None not in outShape
766 else:
767 raise ValueError("Could not determine the shape of a matrix "
768 "given as '{}' because the expected size of {} contains a "
769 "wildcard. Try to give the shape explicitly with the "
770 "string.".format(value, glyphs.shape(shape)))
771 else:
772 outShape = shape
774 # Create the base matrix.
775 if base.startswith("e_"):
776 position = base[2:].split(",")
777 if len(position) == 1:
778 index = int(position[0]) - 1
779 i, j = index % outShape[0], index // outShape[0]
780 else:
781 i, j = int(position[0]) - 1, int(position[1]) - 1
783 if i >= outShape[0] or j >= outShape[1]:
784 raise ValueError("Out-of-boundary unit at row {}, column {} "
785 "in matrix of shape {} given as '{}'."
786 .format(i + 1, j + 1, glyphs.shape(outShape), value))
788 value = load_sparse([1.0], [i], [j], outShape, typecode)
790 if outShape[1] == 1:
791 assert j == 0
792 string = "e_{}".format(i + 1)
793 elif outShape[0] == 1:
794 assert i == 0
795 string = glyphs.transp("e_{}".format(j + 1))
796 else:
797 string = "e_{},{}".format(i + 1, j + 1)
798 elif base.startswith("|"):
799 element = scalar(base[1:-1])
801 # Pull the factor inside the matrix.
802 element *= factor
803 factor = 1
805 value = load_dense(element, outShape, typecode)
806 string = glyphs.scalar(element)
807 elif base == "I":
808 if outShape[0] != outShape[1]:
809 raise ValueError("Cannot create a non-square identy matrix.")
811 n = outShape[0]
813 value = load_sparse([1.0]*n, range(n), range(n), outShape, typecode)
814 string = glyphs.scalar(1) if n == 1 else glyphs.idmatrix()
815 else:
816 assert False, "Unexpected matrix base string '{}'.".format(base)
818 # Apply a coefficient.
819 if factor != 1:
820 value = value * factor
821 string = glyphs.mul(glyphs.scalar(factor), string)
823 # Finalize the string.
824 if base.startswith("e_"):
825 string = glyphs.matrix(
826 glyphs.compsep(glyphs.shape(outShape), string))
827 elif base.startswith("|"):
828 if outShape != (1, 1):
829 string = glyphs.matrix(string)
830 elif base == "I":
831 pass
832 else:
833 assert False, "Unexpected matrix base string."
835 # Convert between dense and sparse representation.
836 if sparse is True and isinstance(value, cvxopt.matrix):
837 value = cvxopt.sparse(value)
838 elif sparse is False and isinstance(value, cvxopt.spmatrix):
839 value = cvxopt.matrix(value)
840 elif isinstance(value, (tuple, list)):
841 if not value:
842 raise ValueError("Cannot parse an empty tuple or list as a matrix.")
844 rows = len(value)
845 cols = None
846 nested = isinstance(value[0], (tuple, list))
848 if nested:
849 # Both outer and inner container must be lists. Unconditionally make
850 # a copy of the outer container so we can replace any inner tuples.
851 value = list(value)
853 for rowNum, row in enumerate(value):
854 if not isinstance(row, (tuple, list)):
855 raise TypeError("Expected a tuple or list for a matrix row "
856 "but found an element of type {}.".format(
857 type(row).__name__))
859 # Make sure row length is consistent.
860 if cols is None:
861 cols = len(row)
863 if not cols:
864 raise TypeError("Cannot parse an empty tuple or list as"
865 " a matrix row.")
866 elif len(row) != cols:
867 raise TypeError("Rows of differing size in a matrix given "
868 "as a tuple or list: {}, {}.".format(cols, len(row)))
870 if isinstance(row, tuple):
871 value[rowNum] = list(row)
873 value = load_dense(value, (cols, rows), typecode).T
874 elif rows == 1:
875 outShape = blend_shapes(shape, (1, 1))
876 return load_data(
877 value[0], outShape, typecode, sparse, alwaysCopy, legacy)
878 else:
879 outShape = (1, rows) if simple_vector_as_row(shape) else (rows, 1)
880 value = load_dense(value, outShape, typecode)
882 # Recurse for further transformations (broadcasting, sparsity).
883 return load_data(value, shape, typecode, sparse, False, legacy)
884 else:
885 raise TypeError("PICOS can't load an object of type {} as a matrix: {}."
886 .format(type(value).__name__, repr(value)))
888 # HACK: Work around CVXOPT not supporting integer sparse matrices: If
889 # integer is requested and the matrix is currently sparse, turn dense.
890 # Note that users may not request both sparsity and integrality.
891 if typecode == "i" and value.typecode == "d":
892 assert not sparse
893 assert isinstance(value, cvxopt.spmatrix)
894 value = load_dense(value, value.size, typecode)
896 if legacy:
897 # Handle the case of broadcasting a scalar for matrix multiplication.
898 assert inShape is not None, "Conversions must define 'inSize'."
899 if inShape == (1, 1) and None in shape and shape != (None, None) \
900 and 1 not in shape:
901 assert 1 in value.size
903 value = cvxopt.spdiag(value)
905 if sparse is False:
906 value = cvxopt.dense(value)
908 scalarString = glyphs.scalar(value[0])
909 string = glyphs.mul(scalarString, glyphs.idmatrix()) \
910 if scalarString != "1" else glyphs.idmatrix()
912 # Validate the output shape and type.
913 assert value.size == blend_shapes(shape, value.size)
914 assert not typecode or value.typecode == typecode
915 if sparse is None:
916 assert isinstance(value, (cvxopt.matrix, cvxopt.spmatrix))
917 elif sparse:
918 assert isinstance(value, cvxopt.spmatrix)
919 else:
920 assert isinstance(value, cvxopt.matrix)
922 # Fallback to a generic matrix string if no better string was found.
923 if string is None:
924 if value.size == (1, 1):
925 string = glyphs.scalar(value[0])
926 else:
927 string = glyphs.matrix(glyphs.shape(value.size))
929 return value, string
932def load_sparse_data(value, shape=None, typecode=None, alwaysCopy=True):
933 """See :func:`~.data.load_data` with ``sparse = True``."""
934 return load_data(value=value, shape=shape, typecode=typecode, sparse=True,
935 alwaysCopy=alwaysCopy)
938def load_dense_data(value, shape=None, typecode=None, alwaysCopy=True):
939 """See :func:`~.data.load_data` with ``sparse = False``."""
940 return load_data(value=value, shape=shape, typecode=typecode, sparse=False,
941 alwaysCopy=alwaysCopy)
944def convert_and_refine_arguments(*which, refine=True, allowNone=False):
945 """Convert selected function arguments to PICOS expressions.
947 If the selected arguments are already PICOS expressions, they are refined
948 unless disabled. If they are not already PICOS expressions, an attempt is
949 made to load them as constant expressions.
951 :Decorator guarantee:
953 All specified arguments are refined PICOS expressions when the function is
954 exectued.
956 :param bool refine:
957 Whether to refine arguments that are already PICOS expressions.
958 :param bool allowNone:
959 Whether :obj:`None` is passed through to the function.
960 """
961 def decorator(func):
962 @functools.wraps(func)
963 def wrapper(*args, **kwargs):
964 from .exp_affine import Constant
965 from .expression import Expression
967 def name():
968 if hasattr(func, "__qualname__"):
969 return func.__qualname__
970 else:
971 return func.__name__
973 callargs = inspect.getcallargs(func, *args, **kwargs)
975 newargs = {}
976 for key, arg in callargs.items():
977 if key not in which:
978 pass
979 elif allowNone and arg is None:
980 pass
981 elif isinstance(arg, Expression):
982 if refine:
983 arg = arg.refined
984 else:
985 try:
986 arg = Constant(arg)
987 except Exception as error:
988 raise TypeError(
989 "Failed to convert argument '{}' of {} to a PICOS "
990 "constant.".format(key, name())) from error
992 newargs[key] = arg
994 return func(**newargs)
995 return wrapper
996 return decorator
999def convert_operands(
1000 sameShape=False, diagBroadcast=False, scalarLHS=False, scalarRHS=False,
1001 rMatMul=False, lMatMul=False, horiCat=False, vertCat=False,
1002 allowNone=False):
1003 """Convert binary operator operands to PICOS expressions.
1005 A decorator for a binary operator that converts any operand that is not
1006 already a PICOS expression or set into a constant one that fulfills the
1007 given shape requirements, if possible. See :func:`~.data.load_data` for
1008 broadcasting and reshaping rules that apply to raw data.
1010 If both operands are already PICOS expressions and at least one of them is
1011 an affine expression, there is a limited set of broadcasting rules to fix
1012 a detected shape mismatch. If this does not succeed, an exception is raised.
1013 If no operand is affine, the operation is performed even if shapes do not
1014 match. The operation is responsible for dealing with this case.
1016 If either operand is a PICOS set, no broadcasting or reshaping is applied as
1017 set instances have, in general, variable dimensionality. If a set type can
1018 *not* have arbitrary dimensions, then it must validate the element's shape
1019 on its own. In particular, no shape requirement may be given when this
1020 decorator is used on a set method.
1022 :Decorator guarantee:
1024 This decorator guarantees to the binary operator using it that only PICOS
1025 expression or set types will be passed as operands and that any affine
1026 expression already has the proper shape for the operation, based on the
1027 decorator arguments.
1029 :Broadcasting rules for affine expressions:
1031 Currently, only scalar affine expressions are broadcasted to the next
1032 smallest matching shape. This is more limited than the broadcasting behavior
1033 when one of the operands is raw data, but it ensures a symmetric behavior in
1034 case both operands are affine. An exception is the case of diagBroadcast,
1035 where a vector affine expression may be extended to a matrix.
1037 :param bool sameShape: Both operands must have the exact same shape.
1038 :param bool diagBroadcast: Both operands must be square matrices of same
1039 shape. If one operand is a square matrix and the other is a scalar or
1040 vector, the latter is put into the diagonal of a square matrix.
1041 :param bool scalarLHS: The left hand side operand must be scalar.
1042 :param bool scalarRHS: The right hand side operand must be scalar.
1043 :param bool rMatMul: The operation has the shape requirements of normal
1044 matrix multiplication with the second operand on the right side.
1045 :param bool lMatMul: The operation has the shape requirements of reversed
1046 matrix multiplication with the second operand on the left side.
1047 :param bool horiCat: The operation has the shape requirements of horizontal
1048 matrix concatenation.
1049 :param bool vertCat: The operation has the shape requirements of vertical
1050 matrix concatenation.
1051 :param bool allowNone: An operand of :obj:`None` is passed as-is.
1053 :raises TypeError: If matching shapes cannot be produced despite one of the
1054 operands being raw data or an affine expression.
1056 .. note::
1057 Matrix multiplication includes scalar by matrix multiplication, so
1058 either operand may remain a scalar.
1059 """
1060 # Fix a redundancy in allowed argument combinations.
1061 if sameShape and (scalarLHS or scalarRHS):
1062 sameShape = False
1063 scalarLHS = True
1064 scalarRHS = True
1066 # Check arguments; only scalarLHS and scalarRHS may appear together.
1067 anyScalar = scalarLHS or scalarRHS
1068 args = (
1069 sameShape, diagBroadcast, anyScalar, lMatMul, rMatMul, horiCat, vertCat)
1070 selected = len([arg for arg in args if arg])
1071 if selected > 1:
1072 assert False, "Conflicting convert_operands arguments."
1074 def decorator(operator):
1075 @functools.wraps(operator)
1076 def wrapper(lhs, rhs, *args, **kwargs):
1077 def fail(reason="Operand shapes do not match.", error=None):
1078 opName = operator.__qualname__ \
1079 if hasattr(operator, "__qualname__") else operator.__name__
1080 lhsName = lhs.string if hasattr(lhs, "string") else repr(lhs)
1081 rhsName = rhs.string if hasattr(rhs, "string") else repr(rhs)
1082 raise TypeError("Invalid operation {}({}, {}): {}".format(
1083 opName, lhsName, rhsName, reason)) from error
1085 def make_lhs_shape(rhsShape):
1086 if sameShape or diagBroadcast:
1087 return rhsShape
1088 elif horiCat:
1089 return (rhsShape[0], None)
1090 elif vertCat:
1091 return (None, rhsShape[1])
1092 elif rMatMul:
1093 return (None, rhsShape[0])
1094 elif lMatMul:
1095 return (rhsShape[1], None)
1096 elif scalarLHS:
1097 return (1, 1)
1098 else:
1099 return None
1101 def make_rhs_shape(lhsShape):
1102 if sameShape or diagBroadcast:
1103 return lhsShape
1104 elif horiCat:
1105 return (lhsShape[0], None)
1106 elif vertCat:
1107 return (None, lhsShape[1])
1108 elif rMatMul:
1109 return (lhsShape[1], None)
1110 elif lMatMul:
1111 return (None, lhsShape[0])
1112 elif scalarRHS:
1113 return (1, 1)
1114 else:
1115 return None
1117 from .exp_affine import Constant
1118 from .exp_biaffine import BiaffineExpression
1119 from .expression import Expression
1120 from .set import Set
1122 lhsIsExpOrSet = isinstance(lhs, (Expression, Set))
1123 rhsIsExpOrSet = isinstance(rhs, (Expression, Set))
1125 lhsIsSet = lhsIsExpOrSet and isinstance(lhs, Set)
1126 rhsIsSet = rhsIsExpOrSet and isinstance(rhs, Set)
1128 if lhsIsSet:
1129 assert not selected, "convert_operands when used on sets may " \
1130 "not pose any shape requirements."
1132 if lhsIsExpOrSet and rhsIsExpOrSet:
1133 lhsIsAffine = isinstance(lhs, BiaffineExpression)
1134 rhsIsAffine = isinstance(rhs, BiaffineExpression)
1136 # If neither expression is biaffine, it's the operation's job to
1137 # deal with it.
1138 if not lhsIsAffine and not rhsIsAffine:
1139 return operator(lhs, rhs, *args, **kwargs)
1141 # If there are no shape requirements, we are done as both are
1142 # already expressions.
1143 if not selected:
1144 return operator(lhs, rhs, *args, **kwargs)
1146 assert not lhsIsSet # Handled earlier.
1148 # Sets have variable shape, so no adjustment is necessary.
1149 if rhsIsSet:
1150 return operator(lhs, rhs, *args, **kwargs)
1152 lhsShape, rhsShape = lhs.shape, rhs.shape
1154 # Check if already matching size.
1155 if (sameShape or horiCat or vertCat) and lhsShape == rhsShape:
1156 return operator(lhs, rhs, *args, **kwargs)
1158 lm, ln = lhs.shape
1159 rm, rn = rhs.shape
1160 lhsSquare = lm == ln
1161 rhsSquare = rm == rn
1163 # Further check if already matching size.
1164 if (diagBroadcast and lhsShape == rhsShape and lhsSquare) \
1165 or (horiCat and lm == rm) \
1166 or (vertCat and ln == rn) \
1167 or (rMatMul and ln == rm) \
1168 or (lMatMul and rn == lm):
1169 return operator(lhs, rhs, *args, **kwargs)
1171 lhsL, rhsL = len(lhs), len(rhs)
1173 # scalarLHS and scalarRHS are the only two shape requirements
1174 # that may appear together, so handle all combinations here.
1175 if scalarLHS and scalarRHS:
1176 if lhsL == 1 and rhsL == 1:
1177 return operator(lhs, rhs, *args, **kwargs)
1178 elif scalarLHS and lhsL == 1:
1179 return operator(lhs, rhs, *args, **kwargs)
1180 elif scalarRHS and rhsL == 1:
1181 return operator(lhs, rhs, *args, **kwargs)
1183 # Matrix multiplication always accepts scalars.
1184 if (rMatMul or lMatMul) and 1 in (lhsL, rhsL):
1185 return operator(lhs, rhs, *args, **kwargs)
1187 # Broadcast an affine scalar on the left side to match.
1188 if lhsIsAffine and lhsL == 1:
1189 if diagBroadcast and rhsSquare:
1190 lhs = lhs.dupdiag(rm)
1191 else:
1192 lhs = lhs.broadcasted(make_lhs_shape(rhsShape))
1193 return operator(lhs, rhs, *args, **kwargs)
1195 # Broadcast an affine scalar on the right side to match.
1196 if rhsIsAffine and rhsL == 1:
1197 if diagBroadcast and lhsSquare:
1198 rhs = rhs.dupdiag(lm)
1199 else:
1200 rhs = rhs.broadcasted(make_rhs_shape(lhsShape))
1201 return operator(lhs, rhs, *args, **kwargs)
1203 # Diagonally broadcast an affine vector on the left hand side.
1204 if diagBroadcast and lhsIsAffine and rhsSquare \
1205 and 1 in lhsShape and rm in lhsShape:
1206 lhs = lhs.diag
1207 assert lhs.shape == rhs.shape
1208 return operator(lhs, rhs, *args, **kwargs)
1210 # Diagonally broadcast an affine vector on the right hand side.
1211 if diagBroadcast and rhsIsAffine and lhsSquare \
1212 and 1 in rhsShape and lm in rhsShape:
1213 rhs = rhs.diag
1214 assert lhs.shape == rhs.shape
1215 return operator(lhs, rhs, *args, **kwargs)
1217 # At least one of the expressions is affine and we didn't find a
1218 # way to fix a detected shape mismatch. It's our job to error.
1219 fail("The operand shapes of {} and {} do not match.".format(
1220 glyphs.shape(lhsShape), glyphs.shape(rhsShape)))
1221 elif lhsIsExpOrSet:
1222 if allowNone and rhs is None:
1223 return operator(lhs, rhs, *args, **kwargs)
1225 rhsShape = None if lhsIsSet else make_rhs_shape(lhs.shape)
1227 if diagBroadcast and lhs.shape[0] != lhs.shape[1]:
1228 fail("Given that the right hand side operand is not a PICOS"
1229 " expression, the left hand side operand must be square"
1230 " for this operation.")
1232 try:
1233 if diagBroadcast:
1234 try:
1235 rhs = Constant(rhs, shape=(1, 1))
1236 except Exception:
1237 try:
1238 rhs = Constant(rhs, shape=(lhs.shape[0], 1))
1239 except Exception:
1240 rhs = Constant(rhs, shape=rhsShape)
1241 else:
1242 rhs = rhs.diag
1243 else:
1244 rhs = rhs.dupdiag(lhs.shape[0])
1245 elif rMatMul or lMatMul:
1246 if lhs.shape == (1, 1):
1247 # Any shape works.
1248 rhs = Constant(rhs)
1249 else:
1250 try: # Try loading as scalar factor first.
1251 rhs = Constant(rhs, shape=(1, 1))
1252 except Exception:
1253 try: # Try loading as a vector.
1254 if lMatMul:
1255 s = (1, lhs.shape[0])
1256 else:
1257 s = (lhs.shape[1], 1)
1259 rhs = Constant(rhs, shape=s)
1260 except TypeError: # Try loading as a matrix.
1261 rhs = Constant(rhs, shape=rhsShape)
1262 else:
1263 rhs = Constant(rhs, shape=rhsShape)
1264 except Exception as error:
1265 fail("Could not load right hand side as a constant of "
1266 "matching shape.", error)
1267 elif rhsIsExpOrSet:
1268 if allowNone and lhs is None:
1269 return operator(lhs, rhs, *args, **kwargs)
1271 lhsShape = None if rhsIsSet else make_lhs_shape(rhs.shape)
1273 if diagBroadcast and rhs.shape[0] != rhs.shape[1]:
1274 fail("Given that the left hand side operand is not a PICOS "
1275 "expression, the right hand side operand must be square"
1276 " for this operation.")
1278 try:
1279 if diagBroadcast:
1280 try:
1281 lhs = Constant(lhs, shape=(1, 1))
1282 except Exception:
1283 try:
1284 lhs = Constant(lhs, shape=(rhs.shape[0], 1))
1285 except Exception:
1286 lhs = Constant(lhs, shape=lhsShape)
1287 else:
1288 lhs = lhs.diag
1289 else:
1290 lhs = lhs.dupdiag(rhs.shape[0])
1291 if rMatMul or lMatMul:
1292 if rhs.shape == (1, 1):
1293 # Any shape works.
1294 lhs = Constant(lhs)
1295 else:
1296 try: # Try loading as scalar factor first.
1297 lhs = Constant(lhs, shape=(1, 1))
1298 except TypeError:
1299 try: # Try loading as a vector.
1300 if rMatMul:
1301 s = (1, rhs.shape[0])
1302 else:
1303 s = (rhs.shape[1], 1)
1305 lhs = Constant(lhs, shape=s)
1306 except TypeError: # Try loading as a matrix.
1307 lhs = Constant(lhs, shape=lhsShape)
1308 else:
1309 lhs = Constant(lhs, shape=lhsShape)
1310 except Exception as error:
1311 fail("Could not load left hand side as a constant of "
1312 "matching shape.", error)
1313 else:
1314 assert False, "convert_operands is supposed to decorate " \
1315 "expression methods, but neither operand is a PICOS " \
1316 "expression or set."
1318 return operator(lhs, rhs, *args, **kwargs)
1319 return wrapper
1320 return decorator
1323def cvxopt_equals(A, B, absTol=None, relTol=None):
1324 """Whether two CVXOPT (sparse) matrices are numerically equal or close.
1326 For every common entry of ``A`` and ``B``, it is sufficient that one of the
1327 two tolerances, ``absTol`` or ``relTol``, is satisfied.
1329 :param float absTol:
1330 Maximum allowed entrywise absolute difference.
1332 :param float relTol:
1333 Maximum allowed entrywise quotient of absolute difference at the entry
1334 divided by the largest absolute value of any entry in both matrices.
1335 """
1336 if A.size != B.size:
1337 return False
1339 Z = A - B
1341 if not Z:
1342 return True
1344 if not absTol and not relTol:
1345 return False
1347 M = max(abs(Z))
1349 if relTol:
1350 N = max(max(abs(A)), max(abs(B)))
1352 if absTol and relTol:
1353 if M > absTol and M / N > relTol:
1354 return False
1355 elif absTol:
1356 if M > absTol:
1357 return False
1358 else:
1359 if M / N > relTol:
1360 return False
1362 return True
1365def cvxopt_maxdiff(A, B):
1366 """Return the largest absolute difference of two (sparse) CVXOPT matrices.
1368 :raises TypeError: If the matrices are not of the same shape.
1369 """
1370 if A.size != B.size:
1371 raise TypeError("The matrices do not have the same shape.")
1373 # Work around "ValueError: max() arg is an empty sequence" for sparse zero.
1374 if not A and not B:
1375 return 0.0
1377 return max(abs(A - B))
1380def cvxopt_hcat(matrices):
1381 """Concatenate the given CVXOPT (sparse) matrices horizontally.
1383 The resulting matrix is sparse if any input matrix is sparse.
1385 :param list matrices: A list of CVXOPT (sparse) matrices.
1386 """
1387 if not isinstance(matrices, list):
1388 matrices = list(matrices)
1390 sparse = any(isinstance(M, cvxopt.spmatrix) for M in matrices)
1392 matrices = [[matrix] for matrix in matrices]
1394 if sparse:
1395 return cvxopt.sparse(matrices)
1396 else:
1397 return cvxopt.matrix(matrices)
1400def cvxopt_vcat(matrices):
1401 """Concatenate the given CVXOPT (sparse) matrices vertically.
1403 The resulting matrix is sparse if any input matrix is sparse.
1405 :param list matrices: A list of CVXOPT (sparse) matrices.
1406 """
1407 if not isinstance(matrices, list):
1408 matrices = list(matrices)
1410 if any(isinstance(M, cvxopt.spmatrix) for M in matrices):
1411 return cvxopt.sparse(matrices)
1412 else:
1413 return cvxopt.matrix(matrices)
1416def cvxopt_hpsd(matrix):
1417 """Whether the given CVXOPT matrix is hermitian positive semidefinite.
1419 Uses :data:`~picos.settings.RELATIVE_HERMITIANNESS_TOLERANCE` and
1420 :data:`~picos.settings.RELATIVE_SEMIDEFINITENESS_TOLERANCE`.
1422 See also :func:`cvxopt_hpd`.
1424 .. warning::
1426 The semidefiniteness tolerance allows negative, near-zero eigenvalues.
1427 """
1428 if not cvxopt_equals(matrix, matrix.H,
1429 relTol=settings.RELATIVE_HERMITIANNESS_TOLERANCE):
1430 return False
1432 eigenvalues = numpy.linalg.eigvalsh(cvx2np(matrix))
1434 minimum = -(max(list(eigenvalues) + [1.0])
1435 * settings.RELATIVE_SEMIDEFINITENESS_TOLERANCE)
1437 return all(ev >= minimum for ev in eigenvalues)
1440def cvxopt_hpd(matrix):
1441 """Whether the given CVXOPT matrix is hermitian positive definite.
1443 Uses :data:`~picos.settings.RELATIVE_HERMITIANNESS_TOLERANCE`.
1445 See also :func:`cvxopt_hpsd`.
1446 """
1447 if not cvxopt_equals(matrix, matrix.H,
1448 relTol=settings.RELATIVE_HERMITIANNESS_TOLERANCE):
1449 return False
1451 eigenvalues = numpy.linalg.eigvalsh(cvx2np(matrix))
1453 return all(ev > 0 for ev in eigenvalues)
1456def cvxopt_inverse(matrix):
1457 """Return the inverse of the given CVXOPT matrix.
1459 :raises ValueError:
1460 If the matrix is not invertible.
1461 """
1462 matrix_np = cvx2np(matrix)
1464 try:
1465 inverse_np = numpy.linalg.inv(matrix_np)
1466 except numpy.linalg.LinAlgError as error:
1467 raise ValueError("Failed to invert a {} CVXOPT matrix using NumPy."
1468 .format(glyphs.shape(matrix.size))) from error
1470 inverse, _ = load_data(inverse_np, matrix.size)
1472 return inverse
1475def cvxopt_principal_root(matrix):
1476 """Return the principal square root of a symmetric positive semidef. matrix.
1478 Given a real symmetric positive (semi)definite CVXOPT input matrix, returns
1479 its unique positive (semi)definite matrix square root.
1481 .. warning::
1483 Does not validate that the input matrix is symmetric positive
1484 semidefinite and will still return a (useless) matrix if it is not.
1485 """
1486 matrix_np = cvx2np(matrix)
1487 U, s, _ = numpy.linalg.svd(matrix_np, hermitian=True)
1488 root_np = numpy.dot(U*(s**0.5), U.T)
1489 root, _ = load_data(root_np, matrix.size)
1491 return root
1494@lru_cache()
1495def cvxopt_K(m, n, typecode="d"):
1496 """The commutation matrix :math:`K_{(m,n)}` as a CVXOPT sparse matrix."""
1497 d = m*n
1498 V = [1]*d
1499 I = range(d)
1500 J = [(k % n)*m + k // n for k in I]
1501 return cvxopt.spmatrix(V, I, J, (d, d), typecode)
1504def sparse_quadruple(A, reshape=None, preT=False, postT=False):
1505 """Return a sparse representation of the given CVXOPT (sparse) matrix.
1507 :param reshape: If set, then :math:`A` is reshaped on the fly.
1508 :param bool preT: Transpose :math:`A` before reshaping.
1509 :param bool postT: Transpose :math:`A` after reshaping.
1511 :returns: A quadruple of values, row indices, column indices, and shape.
1512 """
1513 if not isinstance(A, (cvxopt.spmatrix, cvxopt.matrix)):
1514 raise TypeError("Input must be a CVXOPT (sparse) matrix.")
1516 m, n = A.size
1518 if reshape:
1519 reshape = load_shape(reshape)
1521 if A.size[0] * A.size[1] != reshape[0] * reshape[1]:
1522 raise TypeError("Cannot reshape from {} to {}.".format(
1523 glyphs.shape(A.size), glyphs.shape(reshape)))
1525 if not A:
1526 V, I, J = [], [], []
1527 p, q = reshape if reshape else (n, m if preT else m, n)
1528 else:
1529 if isinstance(A, cvxopt.matrix):
1530 A = cvxopt.sparse(A)
1532 V, I, J = A.V, A.I, A.J
1534 if preT:
1535 I, J, m, n = J, I, n, m
1537 if reshape:
1538 p, q = reshape
1539 I, J = zip(*[((i + j*m) % p, (i + j*m) // p) for i, j in zip(I, J)])
1540 else:
1541 p, q = m, n
1543 if postT:
1544 I, J, p, q = J, I, q, p
1546 return V, I, J, (p, q)
1549def left_kronecker_I(A, k, reshape=None, preT=False, postT=False):
1550 r"""Return :math:`I_k \otimes A` for a CVXOPT (sparse) matrix :math:`A`.
1552 In other words, if :math:`A` is a :math:`m \times n` CVXOPT (sparse)
1553 matrix, returns a :math:`km \times kn` CVXOPT sparse block matrix with
1554 all blocks of size :math:`m \times n`, the diagonal blocks (horizontal
1555 block index equal to vertical block index) equal to :math:`A`, and all
1556 other blocks zero.
1558 :param reshape: If set, then :math:`A` is reshaped on the fly.
1559 :param bool preT: Transpose :math:`A` before reshaping.
1560 :param bool postT: Transpose :math:`A` after reshaping.
1562 :returns:
1563 If :math:`A` is dense and :math:`k = 1`, a
1564 :func:`CVXOPT dense matrix <cvxopt:cvxopt.matrix>`, otherwise a
1565 :func:`CVXOPT sparse matrix <cvxopt:cvxopt.spmatrix>`.
1566 """
1567 A = A.T if preT else A
1569 if reshape:
1570 A = A if preT else A[:] # Copy if not already fresh.
1571 A.size = reshape
1573 A = A.T if postT else A
1575 if k == 1:
1576 if not preT and not reshape and not postT:
1577 return A + 0 # Always return a copy.
1578 else:
1579 return A
1581 if A.size[0] == A.size[1]:
1582 A = cvxopt.spdiag([A for _ in range(k)])
1583 else:
1584 Z = cvxopt.spmatrix([], [], [], A.size)
1585 A = cvxopt.sparse([[Z]*i + [A] + [Z]*(k-i-1) for i in range(k)])
1587 return A.T if postT else A
1590def right_kronecker_I(A, k, reshape=None, preT=False, postT=False):
1591 r"""Return :math:`A \otimes I_k` for a CVXOPT (sparse) matrix :math:`A`.
1593 :param reshape: If set, then :math:`A` is reshaped on the fly.
1594 :param bool preT: Transpose :math:`A` before reshaping.
1595 :param bool postT: Transpose :math:`A` after reshaping.
1597 :returns:
1598 If :math:`A` is dense and :math:`k = 1`, a
1599 :func:`CVXOPT dense matrix <cvxopt:cvxopt.matrix>`, otherwise a
1600 :func:`CVXOPT sparse matrix <cvxopt:cvxopt.spmatrix>`.
1601 """
1602 if isinstance(A, cvxopt.matrix):
1603 # Dense case: Use NumPy.
1604 A = numpy.array(A)
1606 A = A.T if preT else A
1607 A = A.reshape(reshape, order="F") if reshape else A
1608 A = A.T if postT else A
1610 A = numpy.kron(A, numpy.eye(k))
1612 return load_data(A, sparse=(k > 1))[0]
1613 else:
1614 # Sparse case: Python implementation.
1615 # This is slower than the NumPy approach in general but can handle
1616 # arbitrarily large matrices given that they are sufficiently sparse.
1617 V, I, J, shape = sparse_quadruple(A, reshape, preT, postT)
1618 m, n = shape
1620 if V:
1621 V, I, J = zip(*[(v, k*i + l, k*j + l)
1622 for v, i, j in zip(V, I, J) for l in range(k)])
1624 return cvxopt.spmatrix(V, I, J, (k*m, k*n))
1627def cvx2np(A, reshape=None):
1628 """Convert a CVXOPT (sparse) matrix to a NumPy two-dimensional array.
1630 :param A: The CVXOPT :func:`dense <cvxopt:cvxopt.matrix>` or
1631 :func:`sparse <cvxopt:cvxopt.spmatrix>` matrix to convert.
1632 :param bool reshape: Optional new shape for the converted matrix.
1634 :returns: Converted :class:`NumPy array <numpy.ndarray>`.
1635 """
1636 assert isinstance(A, (cvxopt.matrix, cvxopt.spmatrix))
1638 if isinstance(A, cvxopt.spmatrix):
1639 A = cvxopt.matrix(A)
1641 if reshape:
1642 shape = load_shape(reshape)
1643 return numpy.reshape(A, shape, "F")
1644 else:
1645 return numpy.array(A)
1648def cvx2csc(A):
1649 """Convert a CVXOPT matrix to a SciPy sparse matrix in CSC format."""
1650 import scipy.sparse
1652 assert isinstance(A, (cvxopt.matrix, cvxopt.spmatrix))
1654 if isinstance(A, cvxopt.spmatrix):
1655 csc = tuple(tuple(x) for x in reversed(A.CCS))
1656 return scipy.sparse.csc_matrix(csc, shape=A.size)
1657 else:
1658 return scipy.sparse.csc_matrix(A)
1661def cvx2csr(A):
1662 """Convert a CVXOPT matrix to a SciPy sparse matrix in CSR format."""
1663 import scipy.sparse
1665 assert isinstance(A, (cvxopt.matrix, cvxopt.spmatrix))
1667 if isinstance(A, cvxopt.spmatrix):
1668 csc = tuple(tuple(x) for x in reversed(A.CCS))
1669 csc = scipy.sparse.csc_matrix(csc, shape=A.size)
1670 return csc.tocsr()
1671 else:
1672 return scipy.sparse.csr_matrix(A)
1675def sp2cvx(A):
1676 """Convert a SciPy sparse matrix to a CVXOPT sparse matrix."""
1677 import scipy.sparse
1679 assert isinstance(A, scipy.sparse.spmatrix)
1681 A = A.tocoo()
1683 V = A.data
1684 I = A.row
1685 J = A.col
1687 if issubclass(A.dtype.type, numpy.complexfloating):
1688 tc = "z"
1689 else:
1690 tc = "d"
1692 return cvxopt.spmatrix(V, I, J, A.shape, tc)
1695def make_fraction(p, denominator_limit):
1696 """Convert a float :math:`p` to a limited precision fraction.
1698 :param float p: The float to convert, may be positive or negative infinity.
1699 :param int denominator_limit: The largest allowed denominator.
1701 :returns tuple: A quadruple ``(num, den, pNew, pStr)`` with ``pNew`` the
1702 limited precision version of :math:`p`, ``pStr`` a string representation
1703 of the fraction, and ``num`` and ``den`` the numerator and the
1704 denominator of the fraction, respectively.
1705 """
1706 # LEGACY: Old tools.tracepow allowed this.
1707 if p in ("inf", "-inf"):
1708 p = float(p)
1710 if p in (float("inf"), float("-inf")):
1711 return p, 1, p, glyphs.neg(glyphs.infty) if p < 0 else glyphs.infty
1713 frac = Fraction(p).limit_denominator(denominator_limit)
1714 num = frac.numerator
1715 den = frac.denominator
1716 pNew = float(num) / float(den)
1718 if den == 1:
1719 pStr = glyphs.scalar(num)
1720 else:
1721 pStr = glyphs.clever_div(glyphs.scalar(num), glyphs.scalar(den))
1723 return num, den, pNew, pStr
1726def value(obj, sparse=None, numpy=False):
1727 """Convert (nested) PICOS objects to their current value.
1729 :param obj:
1730 Either a single (PICOS) object that has a ``value`` attribute, such as a
1731 :class:`mutable <picos.expressions.mutable.Mutable>`,
1732 :class:`expression <picos.expressions.expression.Expression>` or
1733 :class:`~picos.modeling.Problem`, or a (nested) :class:`list`,
1734 :class:`tuple` or :class:`dict` thereof.
1736 :param sparse:
1737 If :obj:`None`, retrieved multidimensional values can be returned as
1738 either CVXOPT :func:`sparse <cvxopt:cvxopt.spmatrix>` or
1739 :func:`dense <cvxopt:cvxopt.matrix>` matrices, whichever PICOS stores
1740 internally. If :obj:`True` or :obj:`False`, multidimensional values are
1741 always returned as sparse or dense types, respectively.
1743 :param bool numpy:
1744 If :obj:`True`, retrieved multidimensional values are returned as a
1745 NumPy :class:`~numpy:numpy.ndarray` instead of a CVXOPT type. May not be
1746 set in combination with ``sparse=True``.
1748 :returns:
1749 An object of the same (nested) structure as ``obj``, with every
1750 occurence of any object with a ``value`` attribute replaced by that
1751 attribute's current numeric value. In the case of dictionaries, only the
1752 dictionary values will be converted.
1754 :raises TypeError:
1755 If some object with a ``value`` attribute has a value that cannot be
1756 converted to a matrix by :func:`~.load_data`. This can only happen if
1757 the object in question is not a PICOS object.
1759 :Example:
1761 >>> from picos import RealVariable, value
1762 >>> from pprint import pprint
1763 >>> x = {key: RealVariable(key) for key in ("foo", "bar")}
1764 >>> x["foo"].value = 2
1765 >>> x["bar"].value = 3
1766 >>> pprint(value(x))
1767 {'bar': 3.0, 'foo': 2.0}
1768 """
1769 if sparse and numpy:
1770 raise ValueError("NumPy does not support sparse matrices.")
1772 if isinstance(obj, tuple):
1773 return tuple(value(inner, sparse, numpy) for inner in obj)
1774 elif isinstance(obj, list):
1775 return [value(inner, sparse, numpy) for inner in obj]
1776 elif isinstance(obj, dict):
1777 return {k: value(v, sparse, numpy) for k, v in obj.items()}
1778 else:
1779 if hasattr(obj, "value"):
1780 val = obj.value
1782 if isinstance(val, (int, float, complex)):
1783 return val
1785 # PICOS objects always return their value as a CVXOPT matrix type,
1786 # but this function may be used on other objects with a value
1787 # attribute. Try to convert their value to a CVXOPT matrix first.
1788 if not isinstance(val, (cvxopt.matrix, cvxopt.spmatrix)):
1789 if numpy:
1790 load_sparse = False
1791 elif sparse:
1792 load_sparse = True
1793 else:
1794 load_sparse = None
1796 val = load_data(val, sparse=load_sparse)[0]
1797 elif isinstance(val, cvxopt.spmatrix) \
1798 and (sparse is False or numpy):
1799 val = cvxopt.dense(val)
1801 if numpy:
1802 import numpy as the_numpy
1803 assert isinstance(val, cvxopt.matrix)
1804 val = the_numpy.array(val)
1806 return val
1807 else:
1808 return obj
1811# --------------------------------------
1812__all__ = api_end(_API_START, globals())