 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

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

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, 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 is None else int(shape),

76 None if shape is None else int(shape))

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 != shape:

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 if baseShape is None else baseShape,

103 defaultShape if baseShape is None else baseShape)

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)

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))

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

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.

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:

259 (<3x1 matrix, tc='i'>, '[3×1]')

261 (<1x3 matrix, tc='i'>, '[1×3]')

262 >>> A = [[1,2,3],

263 ... [4,5,6]]

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]

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

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 is None and shape is not None and shape > 1

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.

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 * inShape

464 outLength = outShape * outShape

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):

478 value = load_dense(value, outShape, typecode)

479 elif inShape == (outShape, 1):

480 # Copy columns horizontally.

481 value = load_dense(list(value)*outShape, outShape, typecode)

482 elif inShape == (1, outShape):

483 # Copy rows vertically.

485 list(value)*outShape, (outShape, outShape), 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, outShape), typecode).T

508 else:

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 * inShape

533 outLength = outShape * outShape

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):

547 if value != 0:

548 value = load_dense(value, outShape, typecode)

549 else:

550 value = load_sparse([], [], [], outShape, typecode)

551 elif inShape == (outShape, 1):

552 # Copy columns horizontally.

553 V = list(value.V)*outShape

554 I = list(value.I)*outShape

555 J = [j for j in range(outShape) for _ in range(len(value))]

556 value = load_sparse(V, I, J, outShape, typecode)

557 elif inShape == (1, outShape):

558 # Copy rows vertically.

559 V = list(value.V)*outShape

560 I = [i for i in range(outShape) for _ in range(len(value))]

561 J = list(value.J)*outShape

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, len(value)

571 I, J = value.I, value.J

572 I = [I[k] % n + J[k]*n for k in range(nnz)]

573 J = *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, len(value)

583 I, J = value.I, value.J

584 J = [I[k]*m + J[k] % m for k in range(nnz)]

585 I = *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, outShape)

597 value = value.T

598 else:

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.

632 outShape = blend_shapes(shape, inShape)

634 # Define shorthands.

635 inLength = inShape * inShape

636 outLength = outShape * outShape

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):

656 value = numpy.full(outShape, value)

657 elif inShape == (outShape, 1):

658 # Copy columns horizontally.

659 value = value.repeat(outShape, 1)

660 elif inShape == (1, outShape):

661 # Copy rows vertically.

662 value = value.repeat(outShape, 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:

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)

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), int(inShape))

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, shape))

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) - 1

767 i, j = index % outShape, index // outShape

768 else:

769 i, j = int(position) - 1, int(position) - 1

771 if i >= outShape or j >= outShape:

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:

779 assert j == 0

780 string = "e_{}".format(i + 1)

781 elif outShape == 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 != outShape:

797 raise ValueError("Cannot create a non-square identy matrix.")

799 n = outShape

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, (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))

864 value, 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)

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)

913 else:

914 string = glyphs.matrix(glyphs.shape(value.size))

916 return value, string

920 """See :func:~.data.load_data with sparse = True."""

921 return load_data(value=value, shape=shape, typecode=typecode, sparse=True,

922 alwaysCopy=alwaysCopy)

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

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, None)

1071 elif vertCat:

1072 return (None, rhsShape)

1073 elif rMatMul:

1074 return (None, rhsShape)

1075 elif lMatMul:

1076 return (rhsShape, 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, None)

1087 elif vertCat:

1088 return (None, lhsShape)

1089 elif rMatMul:

1090 return (lhsShape, None)

1091 elif lMatMul:

1092 return (None, lhsShape)

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

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:

1168 return operator(lhs, rhs, *args, **kwargs)

1170 # Broadcast an affine scalar on the right side to match.

1171 if rhsIsAffine and rhsL == 1:

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 = *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)

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:

1429 if A.size * A.size != reshape * reshape:

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 == A.size:

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))

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:

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:

1670 elif sparse:

1672 else:

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())