Coverage for picos/expressions/data.py : 68.57%
Hot-keys on this page
r m x p toggle line displays
j k next/prev highlighted chunk
0 (zero) top of page
1 (one) first highlighted chunk
1# ------------------------------------------------------------------------------
2# Copyright (C) 2019-2020 Maximilian Stahlberg
3# Based on the original picos.expressions module by Guillaume Sagnol.
4#
5# This file is part of PICOS.
6#
7# PICOS is free software: you can redistribute it and/or modify it under the
8# terms of the GNU General Public License as published by the Free Software
9# Foundation, either version 3 of the License, or (at your option) any later
10# version.
11#
12# PICOS is distributed in the hope that it will be useful, but WITHOUT ANY
13# WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR
14# A PARTICULAR PURPOSE. See the GNU General Public License for more details.
15#
16# You should have received a copy of the GNU General Public License along with
17# this program. If not, see <http://www.gnu.org/licenses/>.
18# ------------------------------------------------------------------------------
20"""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 load_data(value, shape=None, typecode=None, sparse=None,
129 alwaysCopy=True, legacy=False):
130 r"""Load a constant numeric data value as a CVXOPT (sparse) matrix.
132 As a user, you never need to call this manually, but you should be aware
133 that PICOS uses this function on any raw data you supply as an operand when
134 working with PICOS expressions. For instance, you can just add a NumPy
135 matrix or an algebraic string such as ``"I"`` to such an expression without
136 worrying about any conversion.
138 :Supported data types:
140 - A NumPy matrix: :class:`numpy.ndarray` or :class:`numpy.matrix`.
141 - A CVXOPT matrix: :obj:`cvxopt.matrix` or :obj:`cvxopt.spmatrix`.
142 - A constant PICOS expression: :class:`~.exp_affine.AffineExpression` or
143 :class:`~.exp_affine.ComplexAffineExpression`.
144 - A Python scalar: :class:`int`, :class:`float` or :class:`complex`.
145 - A flat :class:`tuple` or :class:`list` containing scalars or a
146 :class:`range`, all representing a column vector.
147 - A nested :class:`tuple` or :class:`list` containing scalars. The outer
148 level represents rows and the inner level represents the rows'
149 entries. Allows you to define a :math:`2 \times 3` matrix like this:
151 .. code-block:: python
153 A = [[1, 2, 3],
154 [4, 5, 6]]
156 - A verbatim string description, with rows separated by newline and
157 columns separated by whitespace. The same :math:`2 \times 3` matrix as
158 above:
160 .. code-block:: python
162 A = '''1 2 3
163 4 5 6'''
165 - An algebraic string description:
167 .. list-table::
168 :widths: 1 99
170 * - ``"|a|"``
171 - A matrix with all entries equal to :math:`a`.
172 * - ``"|a|(m,n)"``
173 - A :math:`m \times n` matrix with all entries equal to :math:`a`.
174 * - ``"e_i,j(m,n)"``
175 - A :math:`m \times n` matrix with a :math:`1` at :math:`(i,j)`,
176 indexed from :math:`(1,1)`` to :math:`(m,n)``.
177 * - ``"e_i(m,n)"``
178 - A :math:`m \times n` matrix with a single :math:`1` on the
179 :math:`i`-th coordinate, indexed from :math:`1` in column-major
180 order.
181 * - ``"I"``
182 - The identity matrix.
183 * - ``"I(n)"``
184 - The :math:`n \times n` identiy matrix.
185 * - ``"a…"``
186 - The matrix given by ``"…"`` but multiplied by :math:`a`.
188 Different matrix operations such as addition or multiplication have
189 different requirements with respect to the operands' shapes. The shape
190 of any PICOS expression involved will always be maintained. But when an
191 operation involves both a PICOS expression and raw data, then PICOS will
192 try to broadcast or reshape the raw data such that the operation can be
193 performed.
195 :Broadcasting and reshaping rules:
197 - An input vector without a second axis (for instance a non-nested
198 :class:`tuple` or :class:`list` or a :class:`range`) is interpreted as
199 a row vector if the target shape is of the form ``(None, n)`` with
200 :math:`n > 1`, otherwise it is interpreted as a column vector.
201 - If the target shape is :obj:`None` or ``(None, None)``, then the
202 input's shape is maintained.
203 - If the target shape contains :obj:`None` exactly once, that occurence
204 is replaced by the respective dimensionality of the input data shape.
205 - A scalar is copied (broadcasted) to fill the target shape.
206 - A column (row) vector is copied horizontally (vertically) if its
207 length matches the target row (column) count.
208 - Reshaping from matrix to vector: A matrix is vectorized in
209 column-major (row-major) order if the target shape is a column (row)
210 vector whose length matches the number of matrix entries.
211 - Reshaping from vector to matrix: A column (row) vector is treated as
212 the column-major (row-major) vectorization of the target matrix if the
213 former's length matches the number of the latter's entries.
214 - All other combinations of input and target shape raise an exception.
215 - When an algebraic string description specifies no shape, then the
216 shape argument must be supplied. When both the string and the shape
217 argument specify a shape, then they must be consistent (no
218 broadcasting or reshaping is applied in this case).
220 :param shape: The shape of the resulting matrix. If the input data is of
221 another shape, broadcasting or reshaping is used if possible, otherwise
222 an exception is raised. An integer is treated as the length of a column
223 vector. If this is :obj:`None`, then the target's shape is the input's.
224 If only the target number of rows or columns is :obj:`None`, then only
225 that quantity is chosen according to the input.
226 :type shape: int or tuple or list or None
228 :param str typecode: The numeric type of the resulting matrix. Either
229 ``'d'`` (float), ``'z'`` (complex) or ``'i'`` (integer). If the input
230 data is not already of this type, then it will be converted if possible.
231 If this is not possible, then an exception is raised. If this is
232 :obj:`None`, then the output type is chosen based on the input data.
234 :param sparse: If :obj:`True`, a sparse matrix is returned. If :obj:`False`,
235 a dense matrix is returned. If :obj:`None`, it depends on the sparsity
236 pattern of the input data. If the typecode argument is ``'i'``, then
237 a value of :obj:`True` is not allowed and the returned matrix is dense.
238 :type sparse: bool or None
240 :param bool alwaysCopy: If :obj:`True`, then a copy of the input data is
241 returned even if it already equals the output data. If :obj:`False`,
242 the input value can be returned as such if it is already a CVXOPT matrix
243 of the desired shape, typecode, and sparsity.
245 :param bool legacy: Be compatible with the old ``retrieve_matrix`` function.
246 In particular, if the target shape contains :obj:`None` exactly once and
247 the input data is scalar, treat this as a matrix multiplication case and
248 return the scalar times an identity matrix of appropriate size.
250 :returns: A :class:`tuple` whose first entry is the loaded matrix and whose
251 second argument is a short string for representing the data within
252 algebraic expression strings.
254 :Example:
256 >>> from picos.expressions.data import load_data
257 >>> # Data as (nested) list:
258 >>> load_data([1,2,3])
259 (<3x1 matrix, tc='i'>, '[3×1]')
260 >>> load_data([[1,2,3]])
261 (<1x3 matrix, tc='i'>, '[1×3]')
262 >>> A = [[1,2,3],
263 ... [4,5,6]]
264 >>> load_data(A)
265 (<2x3 matrix, tc='i'>, '[2×3]')
266 >>> # Data as string:
267 >>> value, string = load_data('e_14(7,2)')
268 >>> print(string)
269 [7×2:e_7,2]
270 >>> print(value) #doctest: +NORMALIZE_WHITESPACE
271 [ 0 0 ]
272 [ 0 0 ]
273 [ 0 0 ]
274 [ 0 0 ]
275 [ 0 0 ]
276 [ 0 0 ]
277 [ 0 1.00e+00]
278 >>> load_data('5.3I', (2,2))
279 (<2x2 sparse matrix, tc='d', nnz=2>, '5.3·I')
280 """
281 from .exp_affine import ComplexAffineExpression
282 from .expression import Expression
284 def load_sparse(V, I, J, shape, typecode):
285 """Create a CVXOPT sparse matrix."""
286 # HACK: Work around CVXOPT not supporting integer sparse matrices:
287 # Create a real sparse matrix for now (will be converted later).
288 # Note that users may not request both sparsity and integrality.
289 typecode = "d" if typecode == "i" else typecode
291 try:
292 if typecode:
293 return cvxopt.spmatrix(V, I, J, shape, typecode)
294 else:
295 return cvxopt.spmatrix(V, I, J, shape)
296 except TypeError as error:
297 # Attempt to convert complex typed but real valued input to a real
298 # typed output matrix.
299 if typecode == "d":
300 realV = [x.real for x in V]
301 if realV == V:
302 try:
303 return cvxopt.spmatrix(realV, I, J, shape, typecode)
304 except Exception:
305 pass
307 raise TypeError(
308 "Failed to create a CVXOPT sparse matrix of shape {} and type "
309 "'{}'.".format(glyphs.shape(shape), typecode)) from error
311 def load_dense(value, shape, typecode):
312 """Create a CVXOPT dense matrix."""
313 try:
314 if typecode:
315 return cvxopt.matrix(value, shape, typecode)
316 else:
317 return cvxopt.matrix(value, shape)
318 except TypeError as error:
319 # Attempt to convert complex (real) typed but real/integer (integer)
320 # valued input to a real/integer (integer) typed output matrix.
321 if typecode in "id":
322 try:
323 complexValue = cvxopt.matrix(value, shape, "z")
324 except Exception:
325 pass
326 else:
327 if not any(complexValue.imag()):
328 if typecode == "d":
329 return complexValue.real()
330 else:
331 realData = list(complexValue.real())
332 intData = [int(x) for x in realData]
333 if intData == realData:
334 try:
335 return cvxopt.matrix(intData, shape, "i")
336 except Exception:
337 pass
339 raise TypeError(
340 "Failed to create a CVXOPT dense matrix of shape {} and type "
341 "'{}'.".format(glyphs.shape(shape), typecode)) from error
343 def simple_vector_as_row(shape):
344 """Whether a single-axis vector should be a row or column vector."""
345 return shape[0] is None and shape[1] is not None and shape[1] > 1
347 def broadcast_error():
348 """Raise a broadcasting related :class:TypeError."""
349 raise TypeError("Cannot broadcast or reshape from {} to {}{}."
350 .format(glyphs.shape(inShape), glyphs.shape(shape), " read as {}"
351 .format(glyphs.shape(outShape)) if shape != outShape else ""))
353 def scalar(x, typecode=None):
354 """Load a scalar as either a complex, float, or int."""
355 x = complex(x)
357 if typecode == "z":
358 return x # We are content with complex.
359 else:
360 if not x.imag:
361 x = x.real
363 if typecode == "d":
364 return x # We are content with real.
365 else:
366 if int(x) == x:
367 return int(x)
368 elif not typecode:
369 return x # We are not strict and real is the best.
370 elif not typecode:
371 return x # We are not strict and complex is the best.
373 raise TypeError(
374 "Cannot cast {} according to typecode '{}'.".format(x, typecode))
376 # Normalize the shape argument to a two-dimensional tuple of int or None,
377 # where None means any dimensionality.
378 shape = load_shape(shape, wildcards=True)
380 # Validate the typecode argument.
381 if typecode is not None and typecode not in "dzi":
382 raise ValueError("Typecode argument not a valid CVXOPT typecode: {}."
383 .format(typecode))
385 # Validate the sparsity argument.
386 if sparse not in (None, True, False):
387 raise ValueError("Sparsity argument must be True, False or None.")
389 # CVXOPT sparse matrices may not be integer typed.
390 if typecode == "i":
391 if sparse:
392 raise TypeError("Sparse integer matrices are not implemented with "
393 "CVXOPT, which PICOS uses as its matrix backed.")
394 else:
395 sparse = False
397 # Conversions must retrieve the input shape for further processing.
398 inShape = None
400 # Allow conversions to specify their own string description.
401 string = None
403 # Convert from range to list.
404 if isinstance(value, range):
405 value = list(value)
407 # Try to refine a PICOS expression to a constant affine expression.
408 if isinstance(value, Expression):
409 value = value.refined
411 # Convert the data to a CVXOPT (sparse) matrix of proper shape and type.
412 if isinstance(value, (int, float, complex, numpy.number)):
413 # Convert NumPy numbers to Python ones.
414 if isinstance(value, numpy.number):
415 if isinstance(value, numpy.integer):
416 value = int(value)
417 elif isinstance(value, numpy.floating):
418 value = float(value)
419 elif isinstance(value, numpy.complexfloating):
420 value = complex(value)
421 else:
422 assert False, "Unexpected NumPy numeric type {}.".format(
423 type(value).__name__)
425 # CVXOPT is limited by the system's word length.
426 if isinstance(value, int):
427 if value > sys.maxsize or value < -sys.maxsize - 1:
428 raise ValueError("The number {} is too large to be loaded by "
429 "PICOS.".format(value))
431 inShape = (1, 1)
432 outShape = blend_shapes(shape, inShape)
434 string = glyphs.scalar(value)
435 if value and outShape != (1, 1): # Don't use the matrix glyph on zero.
436 string = glyphs.matrix(string)
438 if not value and sparse is not False:
439 value = load_sparse([], [], [], outShape, typecode)
440 else:
441 value = load_dense(value, outShape, typecode)
443 if sparse:
444 value = cvxopt.sparse(value)
445 elif isinstance(value, cvxopt.matrix):
446 # NOTE: Since the input is already a CVXOPT data type, it is possible
447 # that it can be returned without modification. The 'alwaysCopy'
448 # parameter, if set to True, prevents this. As soon as a copy is
449 # made during processing of the input we set 'alwaysCopy' to False
450 # to prevent an unnecessary second copy of the processed input to
451 # be made later on.
453 # If sparse is requested, let CVXOPT handle the transformation first.
454 if sparse:
455 value = cvxopt.sparse(value)
456 return load_data(value, shape, typecode, sparse, False, legacy)
458 # Refine the output shape.
459 inShape = value.size
460 outShape = blend_shapes(shape, inShape)
462 # Define shorthands.
463 inLength = inShape[0] * inShape[1]
464 outLength = outShape[0] * outShape[1]
466 # If the input is stored as complex and no complex output is explicitly
467 # requested, try to remove an imaginary part of zero.
468 if value.typecode == "z" and typecode != "z" and not any(value.imag()):
469 value = value.real()
470 alwaysCopy = False
472 # Broadcast or reshape the data if necessary and possible.
473 reshaped = True
474 if inShape == outShape:
475 reshaped = False
476 elif inShape == (1, 1):
477 # Broadcast a scalar.
478 value = load_dense(value[0], outShape, typecode)
479 elif inShape == (outShape[0], 1):
480 # Copy columns horizontally.
481 value = load_dense(list(value)*outShape[1], outShape, typecode)
482 elif inShape == (1, outShape[1]):
483 # Copy rows vertically.
484 value = load_dense(
485 list(value)*outShape[0], (outShape[1], outShape[0]), typecode).T
486 elif (inLength, 1) == outShape:
487 # Vectorize in column-major order.
488 if not typecode or value.typecode == typecode:
489 # Quick method.
490 value = value[:]
491 else:
492 # With typecode change.
493 value = load_dense(value, outShape, typecode)
494 elif (1, inLength) == outShape:
495 # Vectorize in row-major order.
496 if not typecode or value.typecode == typecode:
497 # Quick method.
498 value = value.T[:].T
499 else:
500 # With typecode change.
501 value = load_dense(value.T, outShape, typecode)
502 elif (outLength, 1) == inShape:
503 # Devectorize in column-major order.
504 value = load_dense(value, outShape, typecode)
505 elif (1, outLength) == inShape:
506 # Devectorize in row-major order.
507 value = load_dense(value, (outShape[1], outShape[0]), typecode).T
508 else:
509 broadcast_error()
510 if reshaped:
511 alwaysCopy = False
513 # The data now has the desired shape.
514 assert value.size == outShape
516 # Ensure the proper typecode.
517 if typecode and value.typecode != typecode:
518 value = load_dense(value, outShape, typecode)
519 alwaysCopy = False
521 # Copy the data if requested and not already a copy.
522 if alwaysCopy:
523 value = load_dense(value, outShape, typecode)
524 elif isinstance(value, cvxopt.spmatrix):
525 # NOTE: See case above.
527 # Refine the output shape.
528 inShape = value.size
529 outShape = blend_shapes(shape, inShape)
531 # Define shorthands.
532 inLength = inShape[0] * inShape[1]
533 outLength = outShape[0] * outShape[1]
535 # If the input is stored as complex and no complex output is explicitly
536 # requested, try to remove an imaginary part of zero.
537 if value.typecode == "z" and typecode != "z" and not any(value.imag()):
538 value = value.real()
539 alwaysCopy = False
541 # Broadcast or reshape the data if necessary and possible.
542 reshaped = True
543 if inShape == outShape:
544 reshaped = False
545 elif inShape == (1, 1):
546 # Broadcast a scalar.
547 if value[0] != 0:
548 value = load_dense(value[0], outShape, typecode)
549 else:
550 value = load_sparse([], [], [], outShape, typecode)
551 elif inShape == (outShape[0], 1):
552 # Copy columns horizontally.
553 V = list(value.V)*outShape[1]
554 I = list(value.I)*outShape[1]
555 J = [j for j in range(outShape[1]) for _ in range(len(value))]
556 value = load_sparse(V, I, J, outShape, typecode)
557 elif inShape == (1, outShape[1]):
558 # Copy rows vertically.
559 V = list(value.V)*outShape[0]
560 I = [i for i in range(outShape[0]) for _ in range(len(value))]
561 J = list(value.J)*outShape[0]
562 value = load_sparse(V, I, J, outShape, typecode)
563 elif (inLength, 1) == outShape:
564 # Vectorize in column-major order.
565 if not typecode or value.typecode == typecode:
566 # Quick method.
567 value = value[:]
568 else:
569 # With typecode change.
570 n, nnz = inShape[0], len(value)
571 I, J = value.I, value.J
572 I = [I[k] % n + J[k]*n for k in range(nnz)]
573 J = [0]*nnz
574 value = load_sparse(value.V, I, J, outShape, typecode)
575 elif (1, inLength) == outShape:
576 # Vectorize in row-major order.
577 if not typecode or value.typecode == typecode:
578 # Quick method.
579 value = value.T[:].T
580 else:
581 # With typecode change.
582 m, nnz = inShape[1], len(value)
583 I, J = value.I, value.J
584 J = [I[k]*m + J[k] % m for k in range(nnz)]
585 I = [0]*nnz
586 value = load_sparse(value.V, I, J, outShape, typecode)
587 elif (outLength, 1) == inShape:
588 # Devectorize in column-major order.
589 # TODO: Logic for also changing the typecode.
590 value = value[:]
591 value.size = outShape
592 elif (1, outLength) == inShape:
593 # Devectorize in row-major order.
594 # TODO: Logic for also changing the typecode.
595 value = value[:]
596 value.size = (outShape[1], outShape[0])
597 value = value.T
598 else:
599 broadcast_error()
600 if reshaped:
601 alwaysCopy = False
603 # The data now has the desired shape.
604 assert value.size == outShape
606 # Ensure the proper typecode.
607 if typecode and value.typecode != typecode:
608 # NOTE: In the case of intential loading as a dense matrix, the
609 # typecode is already set properly.
610 assert isinstance(value, cvxopt.spmatrix)
611 value = load_sparse(value.V, value.I, value.J, outShape, typecode)
612 alwaysCopy = False
614 # Return either a (sparse) copy or the input data itself.
615 if isinstance(value, cvxopt.matrix):
616 # The data was intentionally copied to a dense matrix (through
617 # broadcasting), so only convert it back to sparse if requested.
618 assert reshaped
619 if sparse:
620 value = cvxopt.sparse(value)
621 elif sparse is False:
622 value = load_dense(value, outShape, typecode)
623 elif alwaysCopy:
624 value = load_sparse(value.V, value.I, value.J, outShape, typecode)
625 elif isinstance(value, numpy.ndarray):
626 # NumPy arrays can be tensors, we don't support those.
627 if len(value.shape) > 2:
628 raise TypeError("PICOS does not support tensor data.")
630 # Refine the output shape.
631 inShape = load_shape(value.shape)
632 outShape = blend_shapes(shape, inShape)
634 # Define shorthands.
635 inLength = inShape[0] * inShape[1]
636 outLength = outShape[0] * outShape[1]
638 # If the input is one-dimensional, turn it into a column or row vector.
639 if inShape != value.shape:
640 if simple_vector_as_row(shape):
641 value = value.reshape((1, inLength))
642 else:
643 value = value.reshape(inShape)
645 # If the input is stored as complex and no complex output is explicitly
646 # requested, try to remove an imaginary part of zero.
647 if value.dtype.kind == "c" and typecode != "z" \
648 and not numpy.iscomplex(value).any():
649 value = value.real
651 # Broadcast or reshape the data if necessary and possible.
652 if inShape == outShape:
653 pass
654 elif inShape == (1, 1):
655 # Broadcast a scalar.
656 value = numpy.full(outShape, value)
657 elif inShape == (outShape[0], 1):
658 # Copy columns horizontally.
659 value = value.repeat(outShape[1], 1)
660 elif inShape == (1, outShape[1]):
661 # Copy rows vertically.
662 value = value.repeat(outShape[0], 0)
663 elif (inLength, 1) == outShape:
664 # Vectorize in column-major order.
665 value = value.reshape(outShape, order="F")
666 elif (1, inLength) == outShape:
667 # Vectorize in row-major order.
668 value = value.reshape(outShape, order="C")
669 elif (outLength, 1) == inShape:
670 # Devectorize in column-major order.
671 value = value.reshape(outShape, order="F")
672 elif (1, outLength) == inShape:
673 # Devectorize in row-major order.
674 value = value.reshape(outShape, order="C")
675 else:
676 broadcast_error()
678 # The data now has the desired shape.
679 assert value.shape == outShape
681 # Decide on whether to create a dense or sparse matrix.
682 outSparse = should_be_sparse(outShape, numpy.count_nonzero(value)) \
683 if sparse is None else sparse
685 # Convert to CVXOPT.
686 if outSparse:
687 I, J = value.nonzero()
688 V = value[I, J]
689 value = load_sparse(V, I, J, outShape, typecode)
690 else:
691 value = load_dense(value, outShape, typecode)
692 elif isinstance(value, ComplexAffineExpression):
693 # Must be constant.
694 if not value.constant:
695 raise ValueError("Cannot load the nonconstant expression {} as a "
696 "constant data value.".format(value.string))
698 # Retrieve a copy of the numeric value.
699 value = value.value
701 # NOTE: alwaysCopy=False as the value is already a copy.
702 return load_data(value, shape, typecode, sparse, False, legacy)
703 elif isinstance(value, str) and re.search(r"\s", value):
704 value = value.strip().splitlines()
705 value = [row.strip().split() for row in value]
706 value = [[scalar(x, typecode) for x in row] for row in value]
708 return load_data(value, shape, typecode, sparse, alwaysCopy, legacy)
709 elif isinstance(value, str):
710 # Verbatim matrix strings with only one element fall throught to this
711 # case as they don't (have to) contain a whitespace character.
712 try:
713 value = scalar(value, typecode)
714 except ValueError:
715 pass
716 else:
717 return load_data(value, shape, typecode, sparse, alwaysCopy, legacy)
719 match = _LOAD_DATA_REGEX.match(value)
721 if not match:
722 raise ValueError("The string '{}' could not be parsed as a matrix."
723 .format(value))
725 # Retrieve the string tokens.
726 tokens = match.groups()
727 assert len(tokens) == 3
728 factor, base, inShape = tokens
729 assert base is not None
731 # Convert the factor.
732 factor = scalar(factor) if factor else 1
734 # Determine whether the matrix will be square.
735 square = base == "I"
737 # Convert the shape.
738 if inShape:
739 inShape = inShape[1:-1].split(",")
740 if len(inShape) == 1:
741 inShape *= 2
742 inShape = (int(inShape[0]), int(inShape[1]))
744 if blend_shapes(shape, inShape) != inShape:
745 raise ValueError(
746 "Inconsistent shapes for matrix given as '{}' with expected"
747 " shape {}.".format(value, glyphs.shape(shape)))
749 outShape = inShape
750 elif None in shape:
751 if square and shape != (None, None):
752 outShape = blend_shapes(shape, (shape[1], shape[0]))
753 assert None not in outShape
754 else:
755 raise ValueError("Could not determine the shape of a matrix "
756 "given as '{}' because the expected size of {} contains a "
757 "wildcard. Try to give the shape explicitly with the "
758 "string.".format(value, glyphs.shape(shape)))
759 else:
760 outShape = shape
762 # Create the base matrix.
763 if base.startswith("e_"):
764 position = base[2:].split(",")
765 if len(position) == 1:
766 index = int(position[0]) - 1
767 i, j = index % outShape[0], index // outShape[0]
768 else:
769 i, j = int(position[0]) - 1, int(position[1]) - 1
771 if i >= outShape[0] or j >= outShape[1]:
772 raise ValueError("Out-of-boundary unit at row {}, column {} "
773 "in matrix of shape {} given as '{}'."
774 .format(i + 1, j + 1, glyphs.shape(outShape), value))
776 value = load_sparse([1.0], [i], [j], outShape, typecode)
778 if outShape[1] == 1:
779 assert j == 0
780 string = "e_{}".format(i + 1)
781 elif outShape[0] == 1:
782 assert i == 0
783 string = glyphs.transp("e_{}".format(j + 1))
784 else:
785 string = "e_{},{}".format(i + 1, j + 1)
786 elif base.startswith("|"):
787 element = scalar(base[1:-1])
789 # Pull the factor inside the matrix.
790 element *= factor
791 factor = 1
793 value = load_dense(element, outShape, typecode)
794 string = glyphs.scalar(element)
795 elif base == "I":
796 if outShape[0] != outShape[1]:
797 raise ValueError("Cannot create a non-square identy matrix.")
799 n = outShape[0]
801 value = load_sparse([1.0]*n, range(n), range(n), outShape, typecode)
802 string = glyphs.scalar(1) if n == 1 else glyphs.idmatrix()
803 else:
804 assert False, "Unexpected matrix base string '{}'.".format(base)
806 # Apply a coefficient.
807 if factor != 1:
808 value = value * factor
809 string = glyphs.mul(glyphs.scalar(factor), string)
811 # Finalize the string.
812 if base.startswith("e_"):
813 string = glyphs.matrix(
814 glyphs.compsep(glyphs.shape(outShape), string))
815 elif base.startswith("|"):
816 string = glyphs.matrix(string)
817 elif base == "I":
818 pass
819 else:
820 assert False, "Unexpected matrix base string."
822 # Convert between dense and sparse representation.
823 if sparse is True and isinstance(value, cvxopt.matrix):
824 value = cvxopt.sparse(value)
825 elif sparse is False and isinstance(value, cvxopt.spmatrix):
826 value = cvxopt.matrix(value)
827 elif isinstance(value, (tuple, list)):
828 if not value:
829 raise ValueError("Cannot parse an empty tuple or list as a matrix.")
831 rows = len(value)
832 cols = None
833 nested = isinstance(value[0], (tuple, list))
835 if nested:
836 # Both outer and inner container must be lists. Unconditionally make
837 # a copy of the outer container so we can replace any inner tuples.
838 value = list(value)
840 for rowNum, row in enumerate(value):
841 if not isinstance(row, (tuple, list)):
842 raise TypeError("Expected a tuple or list for a matrix row "
843 "but found an element of type {}.".format(
844 type(row).__name__))
846 # Make sure row length is consistent.
847 if cols is None:
848 cols = len(row)
850 if not cols:
851 raise TypeError("Cannot parse an empty tuple or list as"
852 " a matrix row.")
853 elif len(row) != cols:
854 raise TypeError("Rows of differing size in a matrix given "
855 "as a tuple or list: {}, {}.".format(cols, len(row)))
857 if isinstance(row, tuple):
858 value[rowNum] = list(row)
860 value = load_dense(value, (cols, rows), typecode).T
861 elif rows == 1:
862 outShape = blend_shapes(shape, (1, 1))
863 return load_data(
864 value[0], outShape, typecode, sparse, alwaysCopy, legacy)
865 else:
866 outShape = (1, rows) if simple_vector_as_row(shape) else (rows, 1)
867 value = load_dense(value, outShape, typecode)
869 # Recurse for further transformations (broadcasting, sparsity).
870 return load_data(value, shape, typecode, sparse, False, legacy)
871 else:
872 raise TypeError("PICOS can't load an object of type {} as a matrix: {}."
873 .format(type(value).__name__, repr(value)))
875 # HACK: Work around CVXOPT not supporting integer sparse matrices: If
876 # integer is requested and the matrix is currently sparse, turn dense.
877 # Note that users may not request both sparsity and integrality.
878 if typecode == "i" and value.typecode == "d":
879 assert not sparse
880 assert isinstance(value, cvxopt.spmatrix)
881 value = load_dense(value, value.size, typecode)
883 if legacy:
884 # Handle the case of broadcasting a scalar for matrix multiplication.
885 assert inShape is not None, "Conversions must define 'inSize'."
886 if inShape == (1, 1) and None in shape and shape != (None, None) \
887 and 1 not in shape:
888 assert 1 in value.size
890 value = cvxopt.spdiag(value)
892 if sparse is False:
893 value = cvxopt.dense(value)
895 scalarString = glyphs.scalar(value[0])
896 string = glyphs.mul(scalarString, glyphs.idmatrix()) \
897 if scalarString != "1" else glyphs.idmatrix()
899 # Validate the output shape and type.
900 assert value.size == blend_shapes(shape, value.size)
901 assert not typecode or value.typecode == typecode
902 if sparse is None:
903 assert isinstance(value, (cvxopt.matrix, cvxopt.spmatrix))
904 elif sparse:
905 assert isinstance(value, cvxopt.spmatrix)
906 else:
907 assert isinstance(value, cvxopt.matrix)
909 # Fallback to a generic matrix string if no better string was found.
910 if string is None:
911 if value.size == (1, 1):
912 string = glyphs.scalar(value[0])
913 else:
914 string = glyphs.matrix(glyphs.shape(value.size))
916 return value, string
919def load_sparse_data(value, shape=None, typecode=None, alwaysCopy=True):
920 """See :func:`~.data.load_data` with ``sparse = True``."""
921 return load_data(value=value, shape=shape, typecode=typecode, sparse=True,
922 alwaysCopy=alwaysCopy)
925def load_dense_data(value, shape=None, typecode=None, alwaysCopy=True):
926 """See :func:`~.data.load_data` with ``sparse = False``."""
927 return load_data(value=value, shape=shape, typecode=typecode, sparse=False,
928 alwaysCopy=alwaysCopy)
931def convert_and_refine_arguments(*which, refine=True, allowNone=False):
932 """Convert selected function arguments to PICOS expressions.
934 If the selected arguments are already PICOS expressions, they are refined
935 unless disabled. If they are not already PICOS expressions, an attempt is
936 made to load them as constant expressions.
938 :Decorator guarantee:
940 All specified arguments are refined PICOS expressions when the function is
941 exectued.
943 :param bool refine:
944 Whether to refine arguments that are already PICOS expressions.
945 :param bool allowNone:
946 Whether :obj:`None` is passed through to the function.
947 """
948 def decorator(func):
949 @functools.wraps(func)
950 def wrapper(*args, **kwargs):
951 from .exp_affine import Constant
952 from .expression import Expression
954 def name():
955 if hasattr(func, "__qualname__"):
956 return func.__qualname__
957 else:
958 return func.__name__
960 callargs = inspect.getcallargs(func, *args, **kwargs)
962 newargs = {}
963 for key, arg in callargs.items():
964 if key not in which:
965 pass
966 elif allowNone and arg is None:
967 pass
968 elif isinstance(arg, Expression):
969 if refine:
970 arg = arg.refined
971 else:
972 try:
973 arg = Constant(arg)
974 except Exception as error:
975 raise TypeError(
976 "Failed to convert argument '{}' of {} to a PICOS "
977 "constant.".format(key, name())) from error
979 newargs[key] = arg
981 return func(**newargs)
982 return wrapper
983 return decorator
986def convert_operands(
987 sameShape=False, scalarLHS=False, scalarRHS=False, rMatMul=False,
988 lMatMul=False, horiCat=False, vertCat=False, allowNone=False):
989 """Convert binary operator operands to PICOS expressions.
991 A decorator for a binary operator that converts any operand that is not
992 already a PICOS expression or set into a constant one that fulfills the
993 given shape requirements, if possible. See :func:`~.data.load_data` for
994 broadcasting and reshaping rules that apply to raw data.
996 If both operands are already PICOS expressions and at least one of them is
997 an affine expression, there is a limited set of broadcasting rules to fix
998 a detected shape mismatch. If this does not succeed, an exception is raised.
999 If no operand is affine, the operation is performed even if shapes do not
1000 match. The operation is responsible for dealing with this case.
1002 If either operand is a PICOS set, no broadcasting or reshaping is applied as
1003 set instances have, in general, variable dimensionality. If a set type can
1004 *not* have arbitrary dimensions, then it must validate the element's shape
1005 on its own. In particular, no shape requirement may be given when this
1006 decorator is used on a set method.
1008 :Decorator guarantee:
1010 This decorator guarantees to the binary operator using it that only PICOS
1011 expression or set types will be passed as operands and that any affine
1012 expression already has the proper shape for the operation, based on the
1013 decorator arguments.
1015 :Broadcasting rules for affine expressions:
1017 Currently, only scalar affine expressions are broadcasted to the next
1018 smallest matching shape. This is more limited than the broadcasting behavior
1019 when one of the operands is raw data, but it ensures a symmetric behavior in
1020 case both operands are affine.
1022 :param bool sameShape: Both operands must have the exact same shape.
1023 :param bool scalarLHS: The left hand side operand must be scalar.
1024 :param bool scalarRHS: The right hand side operand must be scalar.
1025 :param bool rMatMul: The operation has the shape requirements of normal
1026 matrix multiplication with the second operand on the right side.
1027 :param bool lMatMul: The operation has the shape requirements of reversed
1028 matrix multiplication with the second operand on the left side.
1029 :param bool horiCat: The operation has the shape requirements of horizontal
1030 matrix concatenation.
1031 :param bool vertCat: The operation has the shape requirements of vertical
1032 matrix concatenation.
1033 :param bool allowNone: An operand of :obj:`None` is passed as-is.
1035 :raises TypeError: If matching shapes cannot be produced despite one of the
1036 operands being raw data or an affine expression.
1038 .. note::
1039 Matrix multiplication includes scalar by matrix multiplication, so
1040 either operand may remain a scalar.
1041 """
1042 # Fix a redundancy in allowed argument combinations.
1043 if sameShape and (scalarLHS or scalarRHS):
1044 sameShape = False
1045 scalarLHS = True
1046 scalarRHS = True
1048 # Check arguments, only scalarLHS and scalarRHS may appear together.
1049 anyScalar = scalarLHS or scalarRHS
1050 selected = len([arg for arg in
1051 (sameShape, anyScalar, lMatMul, rMatMul, horiCat, vertCat) if arg])
1052 if selected > 1:
1053 assert False, "Conflicting convert_operands arguments."
1055 def decorator(operator):
1056 @functools.wraps(operator)
1057 def wrapper(lhs, rhs, *args, **kwargs):
1058 def fail(reason="Operand shapes do not match.", error=None):
1059 opName = operator.__qualname__ \
1060 if hasattr(operator, "__qualname__") else operator.__name__
1061 lhsName = lhs.string if hasattr(lhs, "string") else repr(lhs)
1062 rhsName = rhs.string if hasattr(rhs, "string") else repr(rhs)
1063 raise TypeError("Invalid operation {}({}, {}): {}".format(
1064 opName, lhsName, rhsName, reason)) from error
1066 def make_lhs_shape(rhsShape):
1067 if sameShape:
1068 return rhsShape
1069 elif horiCat:
1070 return (rhsShape[0], None)
1071 elif vertCat:
1072 return (None, rhsShape[1])
1073 elif rMatMul:
1074 return (None, rhsShape[0])
1075 elif lMatMul:
1076 return (rhsShape[1], None)
1077 elif scalarLHS:
1078 return (1, 1)
1079 else:
1080 return None
1082 def make_rhs_shape(lhsShape):
1083 if sameShape:
1084 return lhsShape
1085 elif horiCat:
1086 return (lhsShape[0], None)
1087 elif vertCat:
1088 return (None, lhsShape[1])
1089 elif rMatMul:
1090 return (lhsShape[1], None)
1091 elif lMatMul:
1092 return (None, lhsShape[0])
1093 elif scalarRHS:
1094 return (1, 1)
1095 else:
1096 return None
1098 from .exp_affine import Constant
1099 from .exp_biaffine import BiaffineExpression
1100 from .expression import Expression
1101 from .set import Set
1103 lhsIsExpOrSet = isinstance(lhs, (Expression, Set))
1104 rhsIsExpOrSet = isinstance(rhs, (Expression, Set))
1106 lhsIsSet = lhsIsExpOrSet and isinstance(lhs, Set)
1107 rhsIsSet = rhsIsExpOrSet and isinstance(rhs, Set)
1109 if lhsIsSet:
1110 assert not selected, "convert_operands when used on sets may " \
1111 "not pose any shape requirements."
1113 if lhsIsExpOrSet and rhsIsExpOrSet:
1114 lhsIsAffine = isinstance(lhs, BiaffineExpression)
1115 rhsIsAffine = isinstance(rhs, BiaffineExpression)
1117 # If neither expression is biaffine, it's the operation's job to
1118 # deal with it.
1119 if not lhsIsAffine and not rhsIsAffine:
1120 return operator(lhs, rhs, *args, **kwargs)
1122 # If there are no shape requirements, we are done as both are
1123 # already expressions.
1124 if not selected:
1125 return operator(lhs, rhs, *args, **kwargs)
1127 assert not lhsIsSet # Handled earlier.
1129 # Sets have variable shape, so no adjustment is necessary.
1130 if rhsIsSet:
1131 return operator(lhs, rhs, *args, **kwargs)
1133 lhsShape, rhsShape = lhs.shape, rhs.shape
1135 # Check if already matching size.
1136 if (sameShape or horiCat or vertCat) and lhsShape == rhsShape:
1137 return operator(lhs, rhs, *args, **kwargs)
1139 lm, ln = lhs.shape
1140 rm, rn = rhs.shape
1142 # Further check if already matching size.
1143 if (horiCat and lm == rm) \
1144 or (vertCat and ln == rn) \
1145 or (rMatMul and ln == rm) \
1146 or (lMatMul and rn == lm):
1147 return operator(lhs, rhs, *args, **kwargs)
1149 lhsL, rhsL = len(lhs), len(rhs)
1151 # scalarLHS and scalarRHS are the only two shape requirements
1152 # that may appear together, so handle all combinations here.
1153 if scalarLHS and scalarRHS:
1154 if lhsL == 1 and rhsL == 1:
1155 return operator(lhs, rhs, *args, **kwargs)
1156 elif scalarLHS and lhsL == 1:
1157 return operator(lhs, rhs, *args, **kwargs)
1158 elif scalarRHS and rhsL == 1:
1159 return operator(lhs, rhs, *args, **kwargs)
1161 # Matrix multiplication always accepts scalars.
1162 if (rMatMul or lMatMul) and 1 in (lhsL, rhsL):
1163 return operator(lhs, rhs, *args, **kwargs)
1165 # Broadcast an affine scalar on the left side to match.
1166 if lhsIsAffine and lhsL == 1:
1167 lhs = lhs.broadcasted(make_lhs_shape(rhsShape))
1168 return operator(lhs, rhs, *args, **kwargs)
1170 # Broadcast an affine scalar on the right side to match.
1171 if rhsIsAffine and rhsL == 1:
1172 rhs = rhs.broadcasted(make_rhs_shape(lhsShape))
1173 return operator(lhs, rhs, *args, **kwargs)
1175 # At least one of the expressions is affine and we didn't find a
1176 # way to fix a detected shape mismatch. It's our job to error.
1177 fail("The operand shapes of {} and {} do not match.".format(
1178 glyphs.shape(lhsShape), glyphs.shape(rhsShape)))
1179 elif lhsIsExpOrSet:
1180 if allowNone and rhs is None:
1181 return operator(lhs, rhs, *args, **kwargs)
1183 rhsShape = None if lhsIsSet else make_rhs_shape(lhs.shape)
1185 try:
1186 if rMatMul or lMatMul:
1187 if lhs.shape == (1, 1):
1188 # Any shape works.
1189 rhs = Constant(rhs)
1190 else:
1191 try:
1192 rhs = Constant(rhs, shape=(1, 1))
1193 except Exception:
1194 rhs = Constant(rhs, shape=rhsShape)
1195 else:
1196 rhs = Constant(rhs, shape=rhsShape)
1197 except Exception as error:
1198 fail("Could not load right hand side as a constant of "
1199 "matching shape.", error)
1200 elif rhsIsExpOrSet:
1201 if allowNone and lhs is None:
1202 return operator(lhs, rhs, *args, **kwargs)
1204 lhsShape = None if rhsIsSet else make_lhs_shape(rhs.shape)
1206 try:
1207 if rMatMul or lMatMul:
1208 if rhs.shape == (1, 1):
1209 # Any shape works.
1210 lhs = Constant(lhs)
1211 else:
1212 try:
1213 lhs = Constant(lhs, shape=(1, 1))
1214 except TypeError:
1215 lhs = Constant(lhs, shape=lhsShape)
1216 else:
1217 lhs = Constant(lhs, shape=lhsShape)
1218 except Exception as error:
1219 fail("Could not load left hand side as a constant of "
1220 "matching shape.", error)
1221 else:
1222 assert False, "convert_operands is supposed to decorate " \
1223 "expression methods, but neither operand is a PICOS " \
1224 "expression or set."
1226 return operator(lhs, rhs, *args, **kwargs)
1227 return wrapper
1228 return decorator
1231def cvxopt_equals(A, B, absTol=None, relTol=None):
1232 """Whether two CVXOPT (sparse) matrices are numerically equal or close.
1234 For every common entry of ``A`` and ``B``, it is sufficient that one of the
1235 two tolerances, ``absTol`` or ``relTol``, is satisfied.
1237 :param float absTol:
1238 Maximum allowed entrywise absolute difference.
1240 :param float relTol:
1241 Maximum allowed entrywise quotient of absolute difference at the entry
1242 divided by the largest absolute value of any entry in both matrices.
1243 """
1244 if A.size != B.size:
1245 return False
1247 Z = A - B
1249 if not Z:
1250 return True
1252 if not absTol and not relTol:
1253 return False
1255 M = max(abs(Z))
1257 if relTol:
1258 N = max(max(abs(A)), max(abs(B)))
1260 if absTol and relTol:
1261 if M > absTol and M / N > relTol:
1262 return False
1263 elif absTol:
1264 if M > absTol:
1265 return False
1266 else:
1267 if M / N > relTol:
1268 return False
1270 return True
1273def cvxopt_maxdiff(A, B):
1274 """Return the largest absolute difference of two (sparse) CVXOPT matrices.
1276 :raises TypeError: If the matrices are not of the same shape.
1277 """
1278 if A.size != B.size:
1279 raise TypeError("The matrices do not have the same shape.")
1281 # Work around "ValueError: max() arg is an empty sequence" for sparse zero.
1282 if not A and not B:
1283 return 0.0
1285 return max(abs(A - B))
1288def cvxopt_hcat(matrices):
1289 """Concatenate the given CVXOPT (sparse) matrices horizontally.
1291 The resulting matrix is sparse if any input matrix is sparse.
1293 :param list matrices: A list of CVXOPT (sparse) matrices.
1294 """
1295 if not isinstance(matrices, list):
1296 matrices = list(matrices)
1298 sparse = any(isinstance(M, cvxopt.spmatrix) for M in matrices)
1300 matrices = [[matrix] for matrix in matrices]
1302 if sparse:
1303 return cvxopt.sparse(matrices)
1304 else:
1305 return cvxopt.matrix(matrices)
1308def cvxopt_vcat(matrices):
1309 """Concatenate the given CVXOPT (sparse) matrices vertically.
1311 The resulting matrix is sparse if any input matrix is sparse.
1313 :param list matrices: A list of CVXOPT (sparse) matrices.
1314 """
1315 if not isinstance(matrices, list):
1316 matrices = list(matrices)
1318 if any(isinstance(M, cvxopt.spmatrix) for M in matrices):
1319 return cvxopt.sparse(matrices)
1320 else:
1321 return cvxopt.matrix(matrices)
1324def cvxopt_hpsd(matrix):
1325 """Whether the given CVXOPT matrix is hermitian positive semidefinite.
1327 Uses :data:`~picos.settings.RELATIVE_HERMITIANNESS_TOLERANCE` and
1328 :data:`~picos.settings.RELATIVE_SEMIDEFINITENESS_TOLERANCE`.
1330 See also :func:`cvxopt_hpd`.
1332 .. warning::
1334 The semidefiniteness tolerance allows negative, near-zero eigenvalues.
1335 """
1336 if not cvxopt_equals(matrix, matrix.H,
1337 relTol=settings.RELATIVE_HERMITIANNESS_TOLERANCE):
1338 return False
1340 eigenvalues = numpy.linalg.eigvalsh(cvx2np(matrix))
1342 minimum = -(max(list(eigenvalues) + [1.0])
1343 * settings.RELATIVE_SEMIDEFINITENESS_TOLERANCE)
1345 return all(ev >= minimum for ev in eigenvalues)
1348def cvxopt_hpd(matrix):
1349 """Whether the given CVXOPT matrix is hermitian positive definite.
1351 Uses :data:`~picos.settings.RELATIVE_HERMITIANNESS_TOLERANCE`.
1353 See also :func:`cvxopt_hpsd`.
1354 """
1355 if not cvxopt_equals(matrix, matrix.H,
1356 relTol=settings.RELATIVE_HERMITIANNESS_TOLERANCE):
1357 return False
1359 eigenvalues = numpy.linalg.eigvalsh(cvx2np(matrix))
1361 return all(ev > 0 for ev in eigenvalues)
1364def cvxopt_inverse(matrix):
1365 """Return the inverse of the given CVXOPT matrix.
1367 :raises ValueError:
1368 If the matrix is not invertible.
1369 """
1370 matrix_np = cvx2np(matrix)
1372 try:
1373 inverse_np = numpy.linalg.inv(matrix_np)
1374 except numpy.linalg.LinAlgError as error:
1375 raise ValueError("Failed to invert a {} CVXOPT matrix using NumPy."
1376 .format(glyphs.shape(matrix.size))) from error
1378 inverse, _ = load_data(inverse_np, matrix.size)
1380 return inverse
1383def cvxopt_principal_root(matrix):
1384 """Return the principal square root of a symmetric positive semidef. matrix.
1386 Given a real symmetric positive (semi)definite CVXOPT input matrix, returns
1387 its unique positive (semi)definite matrix square root.
1389 .. warning::
1391 Does not validate that the input matrix is symmetric positive
1392 semidefinite and will still return a (useless) matrix if it is not.
1393 """
1394 matrix_np = numpy.array(cvx2np(matrix)) # TODO: GitLab issue #239.
1395 U, s, _ = numpy.linalg.svd(matrix_np, hermitian=True)
1396 root_np = numpy.dot(U*(s**0.5), U.T)
1397 root, _ = load_data(root_np, matrix.size)
1399 return root
1402@lru_cache()
1403def cvxopt_K(m, n, typecode="d"):
1404 """The commutation matrix :math:`K_{(m,n)}` as a CVXOPT sparse matrix."""
1405 d = m*n
1406 V = [1]*d
1407 I = range(d)
1408 J = [(k % n)*m + k // n for k in I]
1409 return cvxopt.spmatrix(V, I, J, (d, d), typecode)
1412def sparse_quadruple(A, reshape=None, preT=False, postT=False):
1413 """Return a sparse representation of the given CVXOPT (sparse) matrix.
1415 :param reshape: If set, then :math:`A` is reshaped on the fly.
1416 :param bool preT: Transpose :math:`A` before reshaping.
1417 :param bool postT: Transpose :math:`A` after reshaping.
1419 :returns: A quadruple of values, row indices, column indices, and shape.
1420 """
1421 if not isinstance(A, (cvxopt.spmatrix, cvxopt.matrix)):
1422 raise TypeError("Input must be a CVXOPT (sparse) matrix.")
1424 m, n = A.size
1426 if reshape:
1427 reshape = load_shape(reshape)
1429 if A.size[0] * A.size[1] != reshape[0] * reshape[1]:
1430 raise TypeError("Cannot reshape from {} to {}.".format(
1431 glyphs.shape(A.size), glyphs.shape(reshape)))
1433 if not A:
1434 V, I, J = [], [], []
1435 p, q = reshape if reshape else (n, m if preT else m, n)
1436 elif isinstance(A, cvxopt.spmatrix):
1437 V, I, J = A.V, A.I, A.J
1439 if preT:
1440 I, J, m, n = J, I, n, m
1442 if reshape:
1443 p, q = reshape
1444 I, J = zip(*[((i + j*m) % p, (i + j*m) // p) for i, j in zip(I, J)])
1445 else:
1446 p, q = m, n
1448 V, I, J = tuple(V), tuple(I), tuple(J)
1449 else:
1450 if reshape:
1451 p, q = reshape
1452 if preT:
1453 V, I, J = zip(*[
1454 (v, (k // m + (k % m)*n) % p, (k // m + (k % m)*n) // p)
1455 for k, v in enumerate(A) if v])
1456 else:
1457 V, I, J = zip(
1458 *[(v, k % p, k // p) for k, v in enumerate(A) if v])
1459 else:
1460 if preT:
1461 p, q = n, m
1462 V, I, J = zip(
1463 *[(v, k // m, k % m) for k, v in enumerate(A) if v])
1464 else:
1465 p, q = m, n
1466 V, I, J = zip(
1467 *[(v, k % m, k // m) for k, v in enumerate(A) if v])
1469 if postT:
1470 I, J, p, q = J, I, q, p
1472 return V, I, J, (p, q)
1475def left_kronecker_I(A, k, reshape=None, preT=False, postT=False):
1476 r"""Return :math:`I_k \otimes A` for a CVXOPT (sparse) matrix :math:`A`.
1478 In other words, if :math:`A` is a :math:`m \times n` CVXOPT (sparse)
1479 matrix, returns a :math:`km \times kn` CVXOPT sparse block matrix with
1480 all blocks of size :math:`m \times n`, the diagonal blocks (horizontal
1481 block index equal to vertical block index) equal to :math:`A`, and all
1482 other blocks zero.
1484 :param reshape: If set, then :math:`A` is reshaped on the fly.
1485 :param bool preT: Transpose :math:`A` before reshaping.
1486 :param bool postT: Transpose :math:`A` after reshaping.
1488 :returns:
1489 If :math:`A` is dense and :math:`k = 1`, a
1490 :func:`CVXOPT dense matrix <cvxopt:cvxopt.matrix>`, otherwise a
1491 :func:`CVXOPT sparse matrix <cvxopt:cvxopt.spmatrix>`.
1492 """
1493 A = A.T if preT else A
1495 if reshape:
1496 A = A if preT else A[:] # Copy if not already fresh.
1497 A.size = reshape
1499 A = A.T if postT else A
1501 if k == 1:
1502 if not preT and not reshape and not postT:
1503 return A + 0 # Always return a copy.
1504 else:
1505 return A
1507 if A.size[0] == A.size[1]:
1508 A = cvxopt.spdiag([A for _ in range(k)])
1509 else:
1510 Z = cvxopt.spmatrix([], [], [], A.size)
1511 A = cvxopt.sparse([[Z]*i + [A] + [Z]*(k-i-1) for i in range(k)])
1513 return A.T if postT else A
1516def right_kronecker_I(A, k, reshape=None, preT=False, postT=False):
1517 r"""Return :math:`A \otimes I_k` for a CVXOPT (sparse) matrix :math:`A`.
1519 :param reshape: If set, then :math:`A` is reshaped on the fly.
1520 :param bool preT: Transpose :math:`A` before reshaping.
1521 :param bool postT: Transpose :math:`A` after reshaping.
1523 :returns:
1524 If :math:`A` is dense and :math:`k = 1`, a
1525 :func:`CVXOPT dense matrix <cvxopt:cvxopt.matrix>`, otherwise a
1526 :func:`CVXOPT sparse matrix <cvxopt:cvxopt.spmatrix>`.
1527 """
1528 if isinstance(A, cvxopt.matrix):
1529 # Dense case: Use NumPy.
1530 A = numpy.array(A)
1532 A = A.T if preT else A
1533 A = A.reshape(reshape, order="F") if reshape else A
1534 A = A.T if postT else A
1536 A = numpy.kron(A, numpy.eye(k))
1538 return load_data(A, sparse=(k > 1))[0]
1539 else:
1540 # Sparse case: Python implementation.
1541 # This is slower than the NumPy approach in general but can handle
1542 # arbitrarily large matrices given that they are sufficiently sparse.
1543 V, I, J, shape = sparse_quadruple(A, reshape, preT, postT)
1544 m, n = shape
1546 if V:
1547 V, I, J = zip(*[(v, k*i + l, k*j + l)
1548 for v, i, j in zip(V, I, J) for l in range(k)])
1550 return cvxopt.spmatrix(V, I, J, (k*m, k*n))
1553def cvx2np(A, reshape=None):
1554 """Convert a CVXOPT (sparse) matrix to a NumPy matrix.
1556 :param A: The CVXOPT :func:`dense <cvxopt:cvxopt.matrix>` or
1557 :func:`sparse <cvxopt:cvxopt.spmatrix>` matrix to convert.
1558 :param bool reshape: Optional new shape for the converted matrix.
1560 :returns: Converted :class:`NumPy matrix <numpy.matrix>`.
1561 """
1562 assert isinstance(A, (cvxopt.matrix, cvxopt.spmatrix))
1564 if isinstance(A, cvxopt.spmatrix):
1565 A = cvxopt.matrix(A)
1567 if reshape:
1568 shape = load_shape(reshape)
1569 return numpy.reshape(A, shape, "F")
1570 else:
1571 return numpy.matrix(A)
1574def make_fraction(p, denominator_limit):
1575 """Convert a float :math:`p` to a limited precision fraction.
1577 :param float p: The float to convert, may be positive or negative infinity.
1578 :param int denominator_limit: The largest allowed denominator.
1580 :returns tuple: A quadruple ``(num, den, pNew, pStr)`` with ``pNew`` the
1581 limited precision version of :math:`p`, ``pStr`` a string representation
1582 of the fraction, and ``num`` and ``den`` the numerator and the
1583 denominator of the fraction, respectively.
1584 """
1585 # LEGACY: Old tools.tracepow allowed this.
1586 if p in ("inf", "-inf"):
1587 p = float(p)
1589 if p in (float("inf"), float("-inf")):
1590 return p, 1, p, glyphs.neg(glyphs.infty) if p < 0 else glyphs.infty
1592 frac = Fraction(p).limit_denominator(denominator_limit)
1593 num = frac.numerator
1594 den = frac.denominator
1595 pNew = float(num) / float(den)
1597 if den == 1:
1598 pStr = glyphs.scalar(num)
1599 else:
1600 pStr = glyphs.clever_div(glyphs.scalar(num), glyphs.scalar(den))
1602 return num, den, pNew, pStr
1605def value(obj, sparse=None, numpy=False):
1606 """Convert (nested) PICOS objects to their current value.
1608 :param obj:
1609 Either a single (PICOS) object that has a ``value`` attribute, such as a
1610 :class:`mutable <picos.expressions.mutable.Mutable>`,
1611 :class:`expression <picos.expressions.expression.Expression>` or
1612 :class:`~picos.modeling.Problem`, or a (nested) :class:`list`,
1613 :class:`tuple` or :class:`dict` thereof.
1615 :param sparse:
1616 If :obj:`None`, retrieved multidimensional values can be returned as
1617 either CVXOPT :func:`sparse <cvxopt:cvxopt.spmatrix>` or
1618 :func:`dense <cvxopt:cvxopt.matrix>` matrices, whichever PICOS stores
1619 internally. If :obj:`True` or :obj:`False`, multidimensional values are
1620 always returned as sparse or dense types, respectively.
1622 :param bool numpy:
1623 If :obj:`True`, retrieved multidimensional values are returned as a
1624 NumPy :class:`~numpy:numpy.ndarray` instead of a CVXOPT type. May not be
1625 set in combination with ``sparse=True``.
1627 :returns:
1628 An object of the same (nested) structure as ``obj``, with every
1629 occurence of any object with a ``value`` attribute replaced by that
1630 attribute's current numeric value. In the case of dictionaries, only the
1631 dictionary values will be converted.
1633 :raises TypeError:
1634 If some object with a ``value`` attribute has a value that cannot be
1635 converted to a matrix by :func:`~.load_data`. This can only happen if
1636 the object in question is not a PICOS object.
1638 :Example:
1640 >>> from picos import RealVariable, value
1641 >>> from pprint import pprint
1642 >>> x = {key: RealVariable(key) for key in ("foo", "bar")}
1643 >>> x["foo"].value = 2
1644 >>> x["bar"].value = 3
1645 >>> pprint(value(x))
1646 {'bar': 3.0, 'foo': 2.0}
1647 """
1648 if sparse and numpy:
1649 raise ValueError("NumPy does not support sparse matrices.")
1651 if isinstance(obj, tuple):
1652 return tuple(value(inner, sparse, numpy) for inner in obj)
1653 elif isinstance(obj, list):
1654 return [value(inner, sparse, numpy) for inner in obj]
1655 elif isinstance(obj, dict):
1656 return {k: value(v, sparse, numpy) for k, v in obj.items()}
1657 else:
1658 if hasattr(obj, "value"):
1659 val = obj.value
1661 if isinstance(val, (int, float, complex)):
1662 return val
1664 # PICOS objects always return their value as a CVXOPT matrix type,
1665 # but this function may be used on other objects with a value
1666 # attribute. Try to convert their value to a CVXOPT matrix first.
1667 if not isinstance(val, (cvxopt.matrix, cvxopt.spmatrix)):
1668 if numpy:
1669 load_sparse = False
1670 elif sparse:
1671 load_sparse = True
1672 else:
1673 load_sparse = None
1675 val = load_data(val, sparse=load_sparse)[0]
1676 elif isinstance(val, cvxopt.spmatrix) \
1677 and (sparse is False or numpy):
1678 val = cvxopt.dense(val)
1680 if numpy:
1681 import numpy as the_numpy
1682 assert isinstance(val, cvxopt.matrix)
1683 val = the_numpy.array(val)
1685 return val
1686 else:
1687 return obj
1690# --------------------------------------
1691__all__ = api_end(_API_START, globals())