Coverage for picos/expressions/data.py: 70.50%

Shortcuts on this page

r m x   toggle line displays

j k   next/prev highlighted chunk

0   (zero) top of page

1   (one) first highlighted chunk

827 statements  

1# ------------------------------------------------------------------------------ 

2# Copyright (C) 2019-2022 Maximilian Stahlberg 

3# Based on the original picos.expressions module by Guillaume Sagnol. 

4# 

5# This file is part of PICOS. 

6# 

7# PICOS is free software: you can redistribute it and/or modify it under the 

8# terms of the GNU General Public License as published by the Free Software 

9# Foundation, either version 3 of the License, or (at your option) any later 

10# version. 

11# 

12# PICOS is distributed in the hope that it will be useful, but WITHOUT ANY 

13# WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR 

14# A PARTICULAR PURPOSE. See the GNU General Public License for more details. 

15# 

16# You should have received a copy of the GNU General Public License along with 

17# this program. If not, see <http://www.gnu.org/licenses/>. 

18# ------------------------------------------------------------------------------ 

19 

20"""Functions to load and work with raw numeric data.""" 

21 

22import functools 

23import inspect 

24import math 

25import re 

26import sys 

27from fractions import Fraction 

28from functools import lru_cache 

29 

30import cvxopt 

31import numpy 

32 

33from .. import glyphs, settings 

34from ..apidoc import api_end, api_start 

35 

36_API_START = api_start(globals()) 

37# ------------------------------- 

38 

39 

40#: Maximum entrywise absolute deviation allowed for numeric equality checks. 

41TOLERANCE = 1e-6 

42 

43 

44def load_shape(shape, squareMatrix=False, wildcards=False): 

45 """Parse the argument as a matrix shape. 

46 

47 PICOS uses this function whenever you supply a shape parameter to a method. 

48 

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

51 

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

73 

74 shape = ( 

75 None if shape[0] is None else int(shape[0]), 

76 None if shape[1] is None else int(shape[1])) 

77 

78 if not wildcards and None in shape: 

79 raise ValueError("Invalid shape (wildcards not allowed): {}." 

80 .format(glyphs.shape(shape))) 

81 

82 if 0 in shape: 

83 raise ValueError("Invalid shape (zero-dimensional axis): {}." 

84 .format(glyphs.shape(shape))) 

85 

86 if squareMatrix and shape[0] != shape[1]: 

87 raise ValueError("Invalid shape for a square matrix: {}." 

88 .format(glyphs.shape(shape))) 

89 

90 return shape 

91 

92 

93def blend_shapes(baseShape, defaultShape): 

94 """Replace wildcards in one shape with entries of the other. 

95 

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

104 

105 

106def should_be_sparse(shape, numNonZero): 

107 """Decide whether a matrix is considered sparse. 

108 

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) 

115 

116 

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

126 

127 

128def is_scipy_spmat(value): 

129 """Report whether value is a SciPy sparse matrix without importing scipy.""" 

130 return hasattr(value, "__module__") \ 

131 and value.__module__.startswith("scipy.sparse") 

132 

133 

134def load_data(value, shape=None, typecode=None, sparse=None, 

135 alwaysCopy=True, legacy=False): 

136 r"""Load a constant numeric data value as a CVXOPT (sparse) matrix. 

137 

138 As a user, you never need to call this manually, but you should be aware 

139 that PICOS uses this function on any raw data you supply as an operand when 

140 working with PICOS expressions. For instance, you can just add a NumPy 

141 matrix or an algebraic string such as ``"I"`` to such an expression without 

142 worrying about any conversion. 

143 

144 :Supported data types: 

145 

146 - A NumPy matrix: :class:`numpy.ndarray` or :class:`numpy.matrix`. 

147 - A SciPy sparse matrix: All from :mod:`scipy.sparse`. 

148 - A CVXOPT matrix: :obj:`cvxopt.matrix` or :obj:`cvxopt.spmatrix`. 

149 - A constant PICOS expression: :class:`~.exp_affine.AffineExpression` or 

150 :class:`~.exp_affine.ComplexAffineExpression`. 

151 - A Python scalar: :class:`int`, :class:`float` or :class:`complex`. 

152 - A flat :class:`tuple` or :class:`list` containing scalars or a 

153 :class:`range`, all representing a column vector. 

154 - A nested :class:`tuple` or :class:`list` containing scalars. The outer 

155 level represents rows and the inner level represents the rows' 

156 entries. Allows you to define a :math:`2 \times 3` matrix like this: 

157 

158 .. code-block:: python 

159 

160 A = [[1, 2, 3], 

161 [4, 5, 6]] 

162 

163 - A verbatim string description, with rows separated by newline and 

164 columns separated by whitespace. The same :math:`2 \times 3` matrix as 

165 above: 

166 

167 .. code-block:: python 

168 

169 A = '''1 2 3 

170 4 5 6''' 

171 

172 - An algebraic string description: 

173 

174 .. list-table:: 

175 :widths: 1 99 

176 

177 * - ``"|a|"`` 

178 - A matrix with all entries equal to :math:`a`. 

179 * - ``"|a|(m,n)"`` 

180 - A :math:`m \times n` matrix with all entries equal to :math:`a`. 

181 * - ``"e_i,j(m,n)"`` 

182 - A :math:`m \times n` matrix with a :math:`1` at :math:`(i,j)`, 

183 indexed from :math:`(1,1)`` to :math:`(m,n)``. 

184 * - ``"e_i(m,n)"`` 

185 - A :math:`m \times n` matrix with a single :math:`1` on the 

186 :math:`i`-th coordinate, indexed from :math:`1` in column-major 

187 order. 

188 * - ``"I"`` 

189 - The identity matrix. 

190 * - ``"I(n)"`` 

191 - The :math:`n \times n` identiy matrix. 

192 * - ``"a…"`` 

193 - The matrix given by ``"…"`` but multiplied by :math:`a`. 

194 

195 Different matrix operations such as addition or multiplication have 

196 different requirements with respect to the operands' shapes. The shape 

197 of any PICOS expression involved will always be maintained. But when an 

198 operation involves both a PICOS expression and raw data, then PICOS will 

199 try to broadcast or reshape the raw data such that the operation can be 

200 performed. 

201 

202 :Broadcasting and reshaping rules: 

203 

204 - An input vector without a second axis (for instance a non-nested 

205 :class:`tuple` or :class:`list` or a :class:`range`) is interpreted as 

206 a row vector if the target shape is of the form ``(None, n)`` with 

207 :math:`n > 1`, otherwise it is interpreted as a column vector. 

208 - If the target shape is :obj:`None` or ``(None, None)``, then the 

209 input's shape is maintained. 

210 - If the target shape contains :obj:`None` exactly once, that occurence 

211 is replaced by the respective dimensionality of the input data shape. 

212 - A scalar is copied (broadcasted) to fill the target shape. 

213 - A column (row) vector is copied horizontally (vertically) if its 

214 length matches the target row (column) count. 

215 - Reshaping from matrix to vector: A matrix is vectorized in 

216 column-major (row-major) order if the target shape is a column (row) 

217 vector whose length matches the number of matrix entries. 

218 - Reshaping from vector to matrix: A column (row) vector is treated as 

219 the column-major (row-major) vectorization of the target matrix if the 

220 former's length matches the number of the latter's entries. 

221 - All other combinations of input and target shape raise an exception. 

222 - When an algebraic string description specifies no shape, then the 

223 shape argument must be supplied. When both the string and the shape 

224 argument specify a shape, then they must be consistent (no 

225 broadcasting or reshaping is applied in this case). 

226 

227 :param shape: The shape of the resulting matrix. If the input data is of 

228 another shape, broadcasting or reshaping is used if possible, otherwise 

229 an exception is raised. An integer is treated as the length of a column 

230 vector. If this is :obj:`None`, then the target's shape is the input's. 

231 If only the target number of rows or columns is :obj:`None`, then only 

232 that quantity is chosen according to the input. 

233 :type shape: int or tuple or list or None 

234 

235 :param str typecode: The numeric type of the resulting matrix. Either 

236 ``'d'`` (float), ``'z'`` (complex) or ``'i'`` (integer). If the input 

237 data is not already of this type, then it will be converted if possible. 

238 If this is not possible, then an exception is raised. If this is 

239 :obj:`None`, then the output type is chosen based on the input data. 

240 

241 :param sparse: If :obj:`True`, a sparse matrix is returned. If :obj:`False`, 

242 a dense matrix is returned. If :obj:`None`, it depends on the sparsity 

243 pattern of the input data. If the typecode argument is ``'i'``, then 

244 a value of :obj:`True` is not allowed and the returned matrix is dense. 

245 :type sparse: bool or None 

246 

247 :param bool alwaysCopy: If :obj:`True`, then a copy of the input data is 

248 returned even if it already equals the output data. If :obj:`False`, 

249 the input value can be returned as such if it is already a CVXOPT matrix 

250 of the desired shape, typecode, and sparsity. 

251 

252 :param bool legacy: Be compatible with the old ``retrieve_matrix`` function. 

253 In particular, if the target shape contains :obj:`None` exactly once and 

254 the input data is scalar, treat this as a matrix multiplication case and 

255 return the scalar times an identity matrix of appropriate size. 

256 

257 :returns: A :class:`tuple` whose first entry is the loaded matrix and whose 

258 second argument is a short string for representing the data within 

259 algebraic expression strings. 

260 

261 :Example: 

262 

263 >>> from picos.expressions.data import load_data 

264 >>> # Data as (nested) list: 

265 >>> load_data([1,2,3]) 

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

267 >>> load_data([[1,2,3]]) 

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

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

270 ... [4,5,6]] 

271 >>> load_data(A) 

272 (<2x3 matrix, tc='i'>, '[2×3]') 

273 >>> # Data as string: 

274 >>> value, string = load_data('e_14(7,2)') 

275 >>> print(string) 

276 [7×2:e_7,2] 

277 >>> print(value) #doctest: +NORMALIZE_WHITESPACE 

278 [ 0 0 ] 

279 [ 0 0 ] 

280 [ 0 0 ] 

281 [ 0 0 ] 

282 [ 0 0 ] 

283 [ 0 0 ] 

284 [ 0 1.00e+00] 

285 >>> load_data('5.3I', (2,2)) 

286 (<2x2 sparse matrix, tc='d', nnz=2>, '5.3·I') 

287 """ 

288 from .exp_affine import ComplexAffineExpression 

289 from .expression import Expression 

290 

291 def load_sparse(V, I, J, shape, typecode): 

292 """Create a CVXOPT sparse matrix.""" 

293 # HACK: Work around CVXOPT not supporting integer sparse matrices: 

294 # Create a real sparse matrix for now (will be converted later). 

295 # Note that users may not request both sparsity and integrality. 

296 typecode = "d" if typecode == "i" else typecode 

297 

298 try: 

299 if typecode: 

300 return cvxopt.spmatrix(V, I, J, shape, typecode) 

301 else: 

302 return cvxopt.spmatrix(V, I, J, shape) 

303 except TypeError as error: 

304 # Attempt to convert complex typed but real valued input to a real 

305 # typed output matrix. 

306 if typecode == "d": 

307 realV = [x.real for x in V] 

308 if realV == V: 

309 try: 

310 return cvxopt.spmatrix(realV, I, J, shape, typecode) 

311 except Exception: 

312 pass 

313 

314 raise TypeError( 

315 "Failed to create a CVXOPT sparse matrix of shape {} and type " 

316 "'{}'.".format(glyphs.shape(shape), typecode)) from error 

317 

318 def load_dense(value, shape, typecode): 

319 """Create a CVXOPT dense matrix.""" 

320 try: 

321 if typecode: 

322 return cvxopt.matrix(value, shape, typecode) 

323 else: 

324 return cvxopt.matrix(value, shape) 

325 except TypeError as error: 

326 # Attempt to convert complex (real) typed but real/integer (integer) 

327 # valued input to a real/integer (integer) typed output matrix. 

328 if typecode in "id": 

329 try: 

330 complexValue = cvxopt.matrix(value, shape, "z") 

331 except Exception: 

332 pass 

333 else: 

334 if not any(complexValue.imag()): 

335 if typecode == "d": 

336 return complexValue.real() 

337 else: 

338 realData = list(complexValue.real()) 

339 intData = [int(x) for x in realData] 

340 if intData == realData: 

341 try: 

342 return cvxopt.matrix(intData, shape, "i") 

343 except Exception: 

344 pass 

345 

346 raise TypeError( 

347 "Failed to create a CVXOPT dense matrix of shape {} and type " 

348 "'{}'.".format(glyphs.shape(shape), typecode)) from error 

349 

350 def simple_vector_as_row(shape): 

351 """Whether a single-axis vector should be a row or column vector.""" 

352 return shape[0] is None and shape[1] is not None and shape[1] > 1 

353 

354 def broadcast_error(): 

355 """Raise a broadcasting related :class:TypeError.""" 

356 raise TypeError("Cannot broadcast or reshape from {} to {}{}." 

357 .format(glyphs.shape(inShape), glyphs.shape(shape), " read as {}" 

358 .format(glyphs.shape(outShape)) if shape != outShape else "")) 

359 

360 def scalar(x, typecode=None): 

361 """Load a scalar as either a complex, float, or int.""" 

362 x = complex(x) 

363 

364 if typecode == "z": 

365 return x # We are content with complex. 

366 else: 

367 if not x.imag: 

368 x = x.real 

369 

370 if typecode == "d": 

371 return x # We are content with real. 

372 else: 

373 if int(x) == x: 

374 return int(x) 

375 elif not typecode: 

376 return x # We are not strict and real is the best. 

377 elif not typecode: 

378 return x # We are not strict and complex is the best. 

379 

380 raise TypeError( 

381 "Cannot cast {} according to typecode '{}'.".format(x, typecode)) 

382 

383 # Normalize the shape argument to a two-dimensional tuple of int or None, 

384 # where None means any dimensionality. 

385 shape = load_shape(shape, wildcards=True) 

386 

387 # Validate the typecode argument. 

388 if typecode is not None and typecode not in "dzi": 

389 raise ValueError("Typecode argument not a valid CVXOPT typecode: {}." 

390 .format(typecode)) 

391 

392 # Validate the sparsity argument. 

393 if sparse not in (None, True, False): 

394 raise ValueError("Sparsity argument must be True, False or None.") 

395 

396 # CVXOPT sparse matrices may not be integer typed. 

397 if typecode == "i": 

398 if sparse: 

399 raise TypeError("Sparse integer matrices are not implemented with " 

400 "CVXOPT, which PICOS uses as its matrix backed.") 

401 else: 

402 sparse = False 

403 

404 # Conversions must retrieve the input shape for further processing. 

405 inShape = None 

406 

407 # Allow conversions to specify their own string description. 

408 string = None 

409 

410 # Convert from range to list. 

411 if isinstance(value, range): 

412 value = list(value) 

413 

414 # Try to refine a PICOS expression to a constant affine expression. 

415 if isinstance(value, Expression): 

416 value = value.refined 

417 

418 # Convert from a SciPy sparse to a CVXOPT sparse matrix. 

419 if is_scipy_spmat(value): 

420 value = sp2cvx(value) 

421 alwaysCopy = False 

422 

423 # Convert the data to a CVXOPT (sparse) matrix of proper shape and type. 

424 if isinstance(value, (int, float, complex, numpy.number)): 

425 # Convert NumPy numbers to Python ones. 

426 if isinstance(value, numpy.number): 

427 if isinstance(value, numpy.integer): 

428 value = int(value) 

429 elif isinstance(value, numpy.floating): 

430 value = float(value) 

431 elif isinstance(value, numpy.complexfloating): 

432 value = complex(value) 

433 else: 

434 assert False, "Unexpected NumPy numeric type {}.".format( 

435 type(value).__name__) 

436 

437 # CVXOPT is limited by the system's word length. 

438 if isinstance(value, int): 

439 if value > sys.maxsize or value < -sys.maxsize - 1: 

440 raise ValueError("The number {} is too large to be loaded by " 

441 "PICOS.".format(value)) 

442 

443 inShape = (1, 1) 

444 outShape = blend_shapes(shape, inShape) 

445 

446 string = glyphs.scalar(value) 

447 if value and outShape != (1, 1): # Don't use the matrix glyph on zero. 

448 string = glyphs.matrix(string) 

449 

450 if not value and sparse is not False: 

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

452 else: 

453 value = load_dense(value, outShape, typecode) 

454 

455 if sparse: 

456 value = cvxopt.sparse(value) 

457 elif isinstance(value, cvxopt.matrix): 

458 # NOTE: Since the input is already a CVXOPT data type, it is possible 

459 # that it can be returned without modification. The 'alwaysCopy' 

460 # parameter, if set to True, prevents this. As soon as a copy is 

461 # made during processing of the input we set 'alwaysCopy' to False 

462 # to prevent an unnecessary second copy of the processed input to 

463 # be made later on. 

464 

465 # If sparse is requested, let CVXOPT handle the transformation first. 

466 if sparse: 

467 value = cvxopt.sparse(value) 

468 return load_data(value, shape, typecode, sparse, False, legacy) 

469 

470 # Refine the output shape. 

471 inShape = value.size 

472 outShape = blend_shapes(shape, inShape) 

473 

474 # Define shorthands. 

475 inLength = inShape[0] * inShape[1] 

476 outLength = outShape[0] * outShape[1] 

477 

478 # If the input is stored as complex and no complex output is explicitly 

479 # requested, try to remove an imaginary part of zero. 

480 if value.typecode == "z" and typecode != "z" and not any(value.imag()): 

481 value = value.real() 

482 alwaysCopy = False 

483 

484 # Broadcast or reshape the data if necessary and possible. 

485 reshaped = True 

486 if inShape == outShape: 

487 reshaped = False 

488 elif inShape == (1, 1): 

489 # Broadcast a scalar. 

490 value = load_dense(value[0], outShape, typecode) 

491 elif inShape == (outShape[0], 1): 

492 # Copy columns horizontally. 

493 value = load_dense(list(value)*outShape[1], outShape, typecode) 

494 elif inShape == (1, outShape[1]): 

495 # Copy rows vertically. 

496 value = load_dense( 

497 list(value)*outShape[0], (outShape[1], outShape[0]), typecode).T 

498 elif (inLength, 1) == outShape: 

499 # Vectorize in column-major order. 

500 if not typecode or value.typecode == typecode: 

501 # Quick method. 

502 value = value[:] 

503 else: 

504 # With typecode change. 

505 value = load_dense(value, outShape, typecode) 

506 elif (1, inLength) == outShape: 

507 # Vectorize in row-major order. 

508 if not typecode or value.typecode == typecode: 

509 # Quick method. 

510 value = value.T[:].T 

511 else: 

512 # With typecode change. 

513 value = load_dense(value.T, outShape, typecode) 

514 elif (outLength, 1) == inShape: 

515 # Devectorize in column-major order. 

516 value = load_dense(value, outShape, typecode) 

517 elif (1, outLength) == inShape: 

518 # Devectorize in row-major order. 

519 value = load_dense(value, (outShape[1], outShape[0]), typecode).T 

520 else: 

521 broadcast_error() 

522 if reshaped: 

523 alwaysCopy = False 

524 

525 # The data now has the desired shape. 

526 assert value.size == outShape 

527 

528 # Ensure the proper typecode. 

529 if typecode and value.typecode != typecode: 

530 value = load_dense(value, outShape, typecode) 

531 alwaysCopy = False 

532 

533 # Copy the data if requested and not already a copy. 

534 if alwaysCopy: 

535 value = load_dense(value, outShape, typecode) 

536 elif isinstance(value, cvxopt.spmatrix): 

537 # NOTE: See case above. 

538 

539 # Refine the output shape. 

540 inShape = value.size 

541 outShape = blend_shapes(shape, inShape) 

542 

543 # Define shorthands. 

544 inLength = inShape[0] * inShape[1] 

545 outLength = outShape[0] * outShape[1] 

546 

547 # If the input is stored as complex and no complex output is explicitly 

548 # requested, try to remove an imaginary part of zero. 

549 if value.typecode == "z" and typecode != "z" and not any(value.imag()): 

550 value = value.real() 

551 alwaysCopy = False 

552 

553 # Broadcast or reshape the data if necessary and possible. 

554 reshaped = True 

555 if inShape == outShape: 

556 reshaped = False 

557 elif inShape == (1, 1): 

558 # Broadcast a scalar. 

559 if value[0] != 0: 

560 value = load_dense(value[0], outShape, typecode) 

561 else: 

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

563 elif inShape == (outShape[0], 1): 

564 # Copy columns horizontally. 

565 V = list(value.V)*outShape[1] 

566 I = list(value.I)*outShape[1] 

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

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

569 elif inShape == (1, outShape[1]): 

570 # Copy rows vertically. 

571 V = list(value.V)*outShape[0] 

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

573 J = list(value.J)*outShape[0] 

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

575 elif (inLength, 1) == outShape: 

576 # Vectorize in column-major order. 

577 if not typecode or value.typecode == typecode: 

578 # Quick method. 

579 value = value[:] 

580 else: 

581 # With typecode change. 

582 n, nnz = inShape[0], len(value) 

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

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

585 J = [0]*nnz 

586 value = load_sparse(value.V, I, J, outShape, typecode) 

587 elif (1, inLength) == outShape: 

588 # Vectorize in row-major order. 

589 if not typecode or value.typecode == typecode: 

590 # Quick method. 

591 value = value.T[:].T 

592 else: 

593 # With typecode change. 

594 m, nnz = inShape[1], len(value) 

595 I, J = value.I, value.J 

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

597 I = [0]*nnz 

598 value = load_sparse(value.V, I, J, outShape, typecode) 

599 elif (outLength, 1) == inShape: 

600 # Devectorize in column-major order. 

601 # TODO: Logic for also changing the typecode. 

602 value = value[:] 

603 value.size = outShape 

604 elif (1, outLength) == inShape: 

605 # Devectorize in row-major order. 

606 # TODO: Logic for also changing the typecode. 

607 value = value[:] 

608 value.size = (outShape[1], outShape[0]) 

609 value = value.T 

610 else: 

611 broadcast_error() 

612 if reshaped: 

613 alwaysCopy = False 

614 

615 # The data now has the desired shape. 

616 assert value.size == outShape 

617 

618 # Ensure the proper typecode. 

619 if typecode and value.typecode != typecode: 

620 # NOTE: In the case of intential loading as a dense matrix, the 

621 # typecode is already set properly. 

622 assert isinstance(value, cvxopt.spmatrix) 

623 value = load_sparse(value.V, value.I, value.J, outShape, typecode) 

624 alwaysCopy = False 

625 

626 # Return either a (sparse) copy or the input data itself. 

627 if isinstance(value, cvxopt.matrix): 

628 # The data was intentionally copied to a dense matrix (through 

629 # broadcasting), so only convert it back to sparse if requested. 

630 assert reshaped 

631 if sparse: 

632 value = cvxopt.sparse(value) 

633 elif sparse is False: 

634 value = load_dense(value, outShape, typecode) 

635 elif alwaysCopy: 

636 value = load_sparse(value.V, value.I, value.J, outShape, typecode) 

637 elif isinstance(value, numpy.ndarray): 

638 # NumPy arrays can be tensors, we don't support those. 

639 if len(value.shape) > 2: 

640 raise TypeError("PICOS does not support tensor data.") 

641 

642 # Refine the output shape. 

643 inShape = load_shape(value.shape) 

644 outShape = blend_shapes(shape, inShape) 

645 

646 # Define shorthands. 

647 inLength = inShape[0] * inShape[1] 

648 outLength = outShape[0] * outShape[1] 

649 

650 # If the input is one-dimensional, turn it into a column or row vector. 

651 if inShape != value.shape: 

652 if simple_vector_as_row(shape): 

653 value = value.reshape((1, inLength)) 

654 else: 

655 value = value.reshape(inShape) 

656 

657 # If the input is stored as complex and no complex output is explicitly 

658 # requested, try to remove an imaginary part of zero. 

659 if value.dtype.kind == "c" and typecode != "z" \ 

660 and not numpy.iscomplex(value).any(): 

661 value = value.real 

662 

663 # Broadcast or reshape the data if necessary and possible. 

664 if inShape == outShape: 

665 pass 

666 elif inShape == (1, 1): 

667 # Broadcast a scalar. 

668 value = numpy.full(outShape, value) 

669 elif inShape == (outShape[0], 1): 

670 # Copy columns horizontally. 

671 value = value.repeat(outShape[1], 1) 

672 elif inShape == (1, outShape[1]): 

673 # Copy rows vertically. 

674 value = value.repeat(outShape[0], 0) 

675 elif (inLength, 1) == outShape: 

676 # Vectorize in column-major order. 

677 value = value.reshape(outShape, order="F") 

678 elif (1, inLength) == outShape: 

679 # Vectorize in row-major order. 

680 value = value.reshape(outShape, order="C") 

681 elif (outLength, 1) == inShape: 

682 # Devectorize in column-major order. 

683 value = value.reshape(outShape, order="F") 

684 elif (1, outLength) == inShape: 

685 # Devectorize in row-major order. 

686 value = value.reshape(outShape, order="C") 

687 else: 

688 broadcast_error() 

689 

690 # The data now has the desired shape. 

691 assert value.shape == outShape 

692 

693 # Decide on whether to create a dense or sparse matrix. 

694 outSparse = should_be_sparse(outShape, numpy.count_nonzero(value)) \ 

695 if sparse is None else sparse 

696 

697 # Convert to CVXOPT. 

698 if outSparse: 

699 I, J = value.nonzero() 

700 V = value[I, J] 

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

702 else: 

703 value = load_dense(value, outShape, typecode) 

704 elif isinstance(value, ComplexAffineExpression): 

705 # Must be constant. 

706 if not value.constant: 

707 raise ValueError("Cannot load the nonconstant expression {} as a " 

708 "constant data value.".format(value.string)) 

709 

710 # Retrieve a copy of the numeric value. 

711 value = value.value 

712 

713 # NOTE: alwaysCopy=False as the value is already a copy. 

714 return load_data(value, shape, typecode, sparse, False, legacy) 

715 elif isinstance(value, str) and re.search(r"\s", value): 

716 value = value.strip().splitlines() 

717 value = [row.strip().split() for row in value] 

718 value = [[scalar(x, typecode) for x in row] for row in value] 

719 

720 return load_data(value, shape, typecode, sparse, alwaysCopy, legacy) 

721 elif isinstance(value, str): 

722 # Verbatim matrix strings with only one element fall throught to this 

723 # case as they don't (have to) contain a whitespace character. 

724 try: 

725 value = scalar(value, typecode) 

726 except ValueError: 

727 pass 

728 else: 

729 return load_data(value, shape, typecode, sparse, alwaysCopy, legacy) 

730 

731 match = _LOAD_DATA_REGEX.match(value) 

732 

733 if not match: 

734 raise ValueError("The string '{}' could not be parsed as a matrix." 

735 .format(value)) 

736 

737 # Retrieve the string tokens. 

738 tokens = match.groups() 

739 assert len(tokens) == 3 

740 factor, base, inShape = tokens 

741 assert base is not None 

742 

743 # Convert the factor. 

744 factor = scalar(factor) if factor else 1 

745 

746 # Determine whether the matrix will be square. 

747 square = base == "I" 

748 

749 # Convert the shape. 

750 if inShape: 

751 inShape = inShape[1:-1].split(",") 

752 if len(inShape) == 1: 

753 inShape *= 2 

754 inShape = (int(inShape[0]), int(inShape[1])) 

755 

756 if blend_shapes(shape, inShape) != inShape: 

757 raise ValueError( 

758 "Inconsistent shapes for matrix given as '{}' with expected" 

759 " shape {}.".format(value, glyphs.shape(shape))) 

760 

761 outShape = inShape 

762 elif None in shape: 

763 if square and shape != (None, None): 

764 outShape = blend_shapes(shape, (shape[1], shape[0])) 

765 assert None not in outShape 

766 else: 

767 raise ValueError("Could not determine the shape of a matrix " 

768 "given as '{}' because the expected size of {} contains a " 

769 "wildcard. Try to give the shape explicitly with the " 

770 "string.".format(value, glyphs.shape(shape))) 

771 else: 

772 outShape = shape 

773 

774 # Create the base matrix. 

775 if base.startswith("e_"): 

776 position = base[2:].split(",") 

777 if len(position) == 1: 

778 index = int(position[0]) - 1 

779 i, j = index % outShape[0], index // outShape[0] 

780 else: 

781 i, j = int(position[0]) - 1, int(position[1]) - 1 

782 

783 if i >= outShape[0] or j >= outShape[1]: 

784 raise ValueError("Out-of-boundary unit at row {}, column {} " 

785 "in matrix of shape {} given as '{}'." 

786 .format(i + 1, j + 1, glyphs.shape(outShape), value)) 

787 

788 value = load_sparse([1.0], [i], [j], outShape, typecode) 

789 

790 if outShape[1] == 1: 

791 assert j == 0 

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

793 elif outShape[0] == 1: 

794 assert i == 0 

795 string = glyphs.transp("e_{}".format(j + 1)) 

796 else: 

797 string = "e_{},{}".format(i + 1, j + 1) 

798 elif base.startswith("|"): 

799 element = scalar(base[1:-1]) 

800 

801 # Pull the factor inside the matrix. 

802 element *= factor 

803 factor = 1 

804 

805 value = load_dense(element, outShape, typecode) 

806 string = glyphs.scalar(element) 

807 elif base == "I": 

808 if outShape[0] != outShape[1]: 

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

810 

811 n = outShape[0] 

812 

813 value = load_sparse([1.0]*n, range(n), range(n), outShape, typecode) 

814 string = glyphs.scalar(1) if n == 1 else glyphs.idmatrix() 

815 else: 

816 assert False, "Unexpected matrix base string '{}'.".format(base) 

817 

818 # Apply a coefficient. 

819 if factor != 1: 

820 value = value * factor 

821 string = glyphs.mul(glyphs.scalar(factor), string) 

822 

823 # Finalize the string. 

824 if base.startswith("e_"): 

825 string = glyphs.matrix( 

826 glyphs.compsep(glyphs.shape(outShape), string)) 

827 elif base.startswith("|"): 

828 if outShape != (1, 1): 

829 string = glyphs.matrix(string) 

830 elif base == "I": 

831 pass 

832 else: 

833 assert False, "Unexpected matrix base string." 

834 

835 # Convert between dense and sparse representation. 

836 if sparse is True and isinstance(value, cvxopt.matrix): 

837 value = cvxopt.sparse(value) 

838 elif sparse is False and isinstance(value, cvxopt.spmatrix): 

839 value = cvxopt.matrix(value) 

840 elif isinstance(value, (tuple, list)): 

841 if not value: 

842 raise ValueError("Cannot parse an empty tuple or list as a matrix.") 

843 

844 rows = len(value) 

845 cols = None 

846 nested = isinstance(value[0], (tuple, list)) 

847 

848 if nested: 

849 # Both outer and inner container must be lists. Unconditionally make 

850 # a copy of the outer container so we can replace any inner tuples. 

851 value = list(value) 

852 

853 for rowNum, row in enumerate(value): 

854 if not isinstance(row, (tuple, list)): 

855 raise TypeError("Expected a tuple or list for a matrix row " 

856 "but found an element of type {}.".format( 

857 type(row).__name__)) 

858 

859 # Make sure row length is consistent. 

860 if cols is None: 

861 cols = len(row) 

862 

863 if not cols: 

864 raise TypeError("Cannot parse an empty tuple or list as" 

865 " a matrix row.") 

866 elif len(row) != cols: 

867 raise TypeError("Rows of differing size in a matrix given " 

868 "as a tuple or list: {}, {}.".format(cols, len(row))) 

869 

870 if isinstance(row, tuple): 

871 value[rowNum] = list(row) 

872 

873 value = load_dense(value, (cols, rows), typecode).T 

874 elif rows == 1: 

875 outShape = blend_shapes(shape, (1, 1)) 

876 return load_data( 

877 value[0], outShape, typecode, sparse, alwaysCopy, legacy) 

878 else: 

879 outShape = (1, rows) if simple_vector_as_row(shape) else (rows, 1) 

880 value = load_dense(value, outShape, typecode) 

881 

882 # Recurse for further transformations (broadcasting, sparsity). 

883 return load_data(value, shape, typecode, sparse, False, legacy) 

884 else: 

885 raise TypeError("PICOS can't load an object of type {} as a matrix: {}." 

886 .format(type(value).__name__, repr(value))) 

887 

888 # HACK: Work around CVXOPT not supporting integer sparse matrices: If 

889 # integer is requested and the matrix is currently sparse, turn dense. 

890 # Note that users may not request both sparsity and integrality. 

891 if typecode == "i" and value.typecode == "d": 

892 assert not sparse 

893 assert isinstance(value, cvxopt.spmatrix) 

894 value = load_dense(value, value.size, typecode) 

895 

896 if legacy: 

897 # Handle the case of broadcasting a scalar for matrix multiplication. 

898 assert inShape is not None, "Conversions must define 'inSize'." 

899 if inShape == (1, 1) and None in shape and shape != (None, None) \ 

900 and 1 not in shape: 

901 assert 1 in value.size 

902 

903 value = cvxopt.spdiag(value) 

904 

905 if sparse is False: 

906 value = cvxopt.dense(value) 

907 

908 scalarString = glyphs.scalar(value[0]) 

909 string = glyphs.mul(scalarString, glyphs.idmatrix()) \ 

910 if scalarString != "1" else glyphs.idmatrix() 

911 

912 # Validate the output shape and type. 

913 assert value.size == blend_shapes(shape, value.size) 

914 assert not typecode or value.typecode == typecode 

915 if sparse is None: 

916 assert isinstance(value, (cvxopt.matrix, cvxopt.spmatrix)) 

917 elif sparse: 

918 assert isinstance(value, cvxopt.spmatrix) 

919 else: 

920 assert isinstance(value, cvxopt.matrix) 

921 

922 # Fallback to a generic matrix string if no better string was found. 

923 if string is None: 

924 if value.size == (1, 1): 

925 string = glyphs.scalar(value[0]) 

926 else: 

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

928 

929 return value, string 

930 

931 

932def load_sparse_data(value, shape=None, typecode=None, alwaysCopy=True): 

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

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

935 alwaysCopy=alwaysCopy) 

936 

937 

938def load_dense_data(value, shape=None, typecode=None, alwaysCopy=True): 

939 """See :func:`~.data.load_data` with ``sparse = False``.""" 

940 return load_data(value=value, shape=shape, typecode=typecode, sparse=False, 

941 alwaysCopy=alwaysCopy) 

942 

943 

944def convert_and_refine_arguments(*which, refine=True, allowNone=False): 

945 """Convert selected function arguments to PICOS expressions. 

946 

947 If the selected arguments are already PICOS expressions, they are refined 

948 unless disabled. If they are not already PICOS expressions, an attempt is 

949 made to load them as constant expressions. 

950 

951 :Decorator guarantee: 

952 

953 All specified arguments are refined PICOS expressions when the function is 

954 exectued. 

955 

956 :param bool refine: 

957 Whether to refine arguments that are already PICOS expressions. 

958 :param bool allowNone: 

959 Whether :obj:`None` is passed through to the function. 

960 """ 

961 def decorator(func): 

962 @functools.wraps(func) 

963 def wrapper(*args, **kwargs): 

964 from .exp_affine import Constant 

965 from .expression import Expression 

966 

967 def name(): 

968 if hasattr(func, "__qualname__"): 

969 return func.__qualname__ 

970 else: 

971 return func.__name__ 

972 

973 callargs = inspect.getcallargs(func, *args, **kwargs) 

974 

975 newargs = {} 

976 for key, arg in callargs.items(): 

977 if key not in which: 

978 pass 

979 elif allowNone and arg is None: 

980 pass 

981 elif isinstance(arg, Expression): 

982 if refine: 

983 arg = arg.refined 

984 else: 

985 try: 

986 arg = Constant(arg) 

987 except Exception as error: 

988 raise TypeError( 

989 "Failed to convert argument '{}' of {} to a PICOS " 

990 "constant.".format(key, name())) from error 

991 

992 newargs[key] = arg 

993 

994 return func(**newargs) 

995 return wrapper 

996 return decorator 

997 

998 

999def convert_operands( 

1000 sameShape=False, diagBroadcast=False, scalarLHS=False, scalarRHS=False, 

1001 rMatMul=False, lMatMul=False, horiCat=False, vertCat=False, 

1002 allowNone=False): 

1003 """Convert binary operator operands to PICOS expressions. 

1004 

1005 A decorator for a binary operator that converts any operand that is not 

1006 already a PICOS expression or set into a constant one that fulfills the 

1007 given shape requirements, if possible. See :func:`~.data.load_data` for 

1008 broadcasting and reshaping rules that apply to raw data. 

1009 

1010 If both operands are already PICOS expressions and at least one of them is 

1011 an affine expression, there is a limited set of broadcasting rules to fix 

1012 a detected shape mismatch. If this does not succeed, an exception is raised. 

1013 If no operand is affine, the operation is performed even if shapes do not 

1014 match. The operation is responsible for dealing with this case. 

1015 

1016 If either operand is a PICOS set, no broadcasting or reshaping is applied as 

1017 set instances have, in general, variable dimensionality. If a set type can 

1018 *not* have arbitrary dimensions, then it must validate the element's shape 

1019 on its own. In particular, no shape requirement may be given when this 

1020 decorator is used on a set method. 

1021 

1022 :Decorator guarantee: 

1023 

1024 This decorator guarantees to the binary operator using it that only PICOS 

1025 expression or set types will be passed as operands and that any affine 

1026 expression already has the proper shape for the operation, based on the 

1027 decorator arguments. 

1028 

1029 :Broadcasting rules for affine expressions: 

1030 

1031 Currently, only scalar affine expressions are broadcasted to the next 

1032 smallest matching shape. This is more limited than the broadcasting behavior 

1033 when one of the operands is raw data, but it ensures a symmetric behavior in 

1034 case both operands are affine. An exception is the case of diagBroadcast, 

1035 where a vector affine expression may be extended to a matrix. 

1036 

1037 :param bool sameShape: Both operands must have the exact same shape. 

1038 :param bool diagBroadcast: Both operands must be square matrices of same 

1039 shape. If one operand is a square matrix and the other is a scalar or 

1040 vector, the latter is put into the diagonal of a square matrix. 

1041 :param bool scalarLHS: The left hand side operand must be scalar. 

1042 :param bool scalarRHS: The right hand side operand must be scalar. 

1043 :param bool rMatMul: The operation has the shape requirements of normal 

1044 matrix multiplication with the second operand on the right side. 

1045 :param bool lMatMul: The operation has the shape requirements of reversed 

1046 matrix multiplication with the second operand on the left side. 

1047 :param bool horiCat: The operation has the shape requirements of horizontal 

1048 matrix concatenation. 

1049 :param bool vertCat: The operation has the shape requirements of vertical 

1050 matrix concatenation. 

1051 :param bool allowNone: An operand of :obj:`None` is passed as-is. 

1052 

1053 :raises TypeError: If matching shapes cannot be produced despite one of the 

1054 operands being raw data or an affine expression. 

1055 

1056 .. note:: 

1057 Matrix multiplication includes scalar by matrix multiplication, so 

1058 either operand may remain a scalar. 

1059 """ 

1060 # Fix a redundancy in allowed argument combinations. 

1061 if sameShape and (scalarLHS or scalarRHS): 

1062 sameShape = False 

1063 scalarLHS = True 

1064 scalarRHS = True 

1065 

1066 # Check arguments; only scalarLHS and scalarRHS may appear together. 

1067 anyScalar = scalarLHS or scalarRHS 

1068 args = ( 

1069 sameShape, diagBroadcast, anyScalar, lMatMul, rMatMul, horiCat, vertCat) 

1070 selected = len([arg for arg in args if arg]) 

1071 if selected > 1: 

1072 assert False, "Conflicting convert_operands arguments." 

1073 

1074 def decorator(operator): 

1075 @functools.wraps(operator) 

1076 def wrapper(lhs, rhs, *args, **kwargs): 

1077 def fail(reason="Operand shapes do not match.", error=None): 

1078 opName = operator.__qualname__ \ 

1079 if hasattr(operator, "__qualname__") else operator.__name__ 

1080 lhsName = lhs.string if hasattr(lhs, "string") else repr(lhs) 

1081 rhsName = rhs.string if hasattr(rhs, "string") else repr(rhs) 

1082 raise TypeError("Invalid operation {}({}, {}): {}".format( 

1083 opName, lhsName, rhsName, reason)) from error 

1084 

1085 def make_lhs_shape(rhsShape): 

1086 if sameShape or diagBroadcast: 

1087 return rhsShape 

1088 elif horiCat: 

1089 return (rhsShape[0], None) 

1090 elif vertCat: 

1091 return (None, rhsShape[1]) 

1092 elif rMatMul: 

1093 return (None, rhsShape[0]) 

1094 elif lMatMul: 

1095 return (rhsShape[1], None) 

1096 elif scalarLHS: 

1097 return (1, 1) 

1098 else: 

1099 return None 

1100 

1101 def make_rhs_shape(lhsShape): 

1102 if sameShape or diagBroadcast: 

1103 return lhsShape 

1104 elif horiCat: 

1105 return (lhsShape[0], None) 

1106 elif vertCat: 

1107 return (None, lhsShape[1]) 

1108 elif rMatMul: 

1109 return (lhsShape[1], None) 

1110 elif lMatMul: 

1111 return (None, lhsShape[0]) 

1112 elif scalarRHS: 

1113 return (1, 1) 

1114 else: 

1115 return None 

1116 

1117 from .exp_affine import Constant 

1118 from .exp_biaffine import BiaffineExpression 

1119 from .expression import Expression 

1120 from .set import Set 

1121 

1122 lhsIsExpOrSet = isinstance(lhs, (Expression, Set)) 

1123 rhsIsExpOrSet = isinstance(rhs, (Expression, Set)) 

1124 

1125 lhsIsSet = lhsIsExpOrSet and isinstance(lhs, Set) 

1126 rhsIsSet = rhsIsExpOrSet and isinstance(rhs, Set) 

1127 

1128 if lhsIsSet: 

1129 assert not selected, "convert_operands when used on sets may " \ 

1130 "not pose any shape requirements." 

1131 

1132 if lhsIsExpOrSet and rhsIsExpOrSet: 

1133 lhsIsAffine = isinstance(lhs, BiaffineExpression) 

1134 rhsIsAffine = isinstance(rhs, BiaffineExpression) 

1135 

1136 # If neither expression is biaffine, it's the operation's job to 

1137 # deal with it. 

1138 if not lhsIsAffine and not rhsIsAffine: 

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

1140 

1141 # If there are no shape requirements, we are done as both are 

1142 # already expressions. 

1143 if not selected: 

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

1145 

1146 assert not lhsIsSet # Handled earlier. 

1147 

1148 # Sets have variable shape, so no adjustment is necessary. 

1149 if rhsIsSet: 

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

1151 

1152 lhsShape, rhsShape = lhs.shape, rhs.shape 

1153 

1154 # Check if already matching size. 

1155 if (sameShape or horiCat or vertCat) and lhsShape == rhsShape: 

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

1157 

1158 lm, ln = lhs.shape 

1159 rm, rn = rhs.shape 

1160 lhsSquare = lm == ln 

1161 rhsSquare = rm == rn 

1162 

1163 # Further check if already matching size. 

1164 if (diagBroadcast and lhsShape == rhsShape and lhsSquare) \ 

1165 or (horiCat and lm == rm) \ 

1166 or (vertCat and ln == rn) \ 

1167 or (rMatMul and ln == rm) \ 

1168 or (lMatMul and rn == lm): 

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

1170 

1171 lhsL, rhsL = len(lhs), len(rhs) 

1172 

1173 # scalarLHS and scalarRHS are the only two shape requirements 

1174 # that may appear together, so handle all combinations here. 

1175 if scalarLHS and scalarRHS: 

1176 if lhsL == 1 and rhsL == 1: 

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

1178 elif scalarLHS and lhsL == 1: 

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

1180 elif scalarRHS and rhsL == 1: 

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

1182 

1183 # Matrix multiplication always accepts scalars. 

1184 if (rMatMul or lMatMul) and 1 in (lhsL, rhsL): 

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

1186 

1187 # Broadcast an affine scalar on the left side to match. 

1188 if lhsIsAffine and lhsL == 1: 

1189 if diagBroadcast and rhsSquare: 

1190 lhs = lhs.dupdiag(rm) 

1191 else: 

1192 lhs = lhs.broadcasted(make_lhs_shape(rhsShape)) 

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

1194 

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

1196 if rhsIsAffine and rhsL == 1: 

1197 if diagBroadcast and lhsSquare: 

1198 rhs = rhs.dupdiag(lm) 

1199 else: 

1200 rhs = rhs.broadcasted(make_rhs_shape(lhsShape)) 

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

1202 

1203 # Diagonally broadcast an affine vector on the left hand side. 

1204 if diagBroadcast and lhsIsAffine and rhsSquare \ 

1205 and 1 in lhsShape and rm in lhsShape: 

1206 lhs = lhs.diag 

1207 assert lhs.shape == rhs.shape 

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

1209 

1210 # Diagonally broadcast an affine vector on the right hand side. 

1211 if diagBroadcast and rhsIsAffine and lhsSquare \ 

1212 and 1 in rhsShape and lm in rhsShape: 

1213 rhs = rhs.diag 

1214 assert lhs.shape == rhs.shape 

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

1216 

1217 # At least one of the expressions is affine and we didn't find a 

1218 # way to fix a detected shape mismatch. It's our job to error. 

1219 fail("The operand shapes of {} and {} do not match.".format( 

1220 glyphs.shape(lhsShape), glyphs.shape(rhsShape))) 

1221 elif lhsIsExpOrSet: 

1222 if allowNone and rhs is None: 

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

1224 

1225 rhsShape = None if lhsIsSet else make_rhs_shape(lhs.shape) 

1226 

1227 if diagBroadcast and lhs.shape[0] != lhs.shape[1]: 

1228 fail("Given that the right hand side operand is not a PICOS" 

1229 " expression, the left hand side operand must be square" 

1230 " for this operation.") 

1231 

1232 try: 

1233 if diagBroadcast: 

1234 try: 

1235 rhs = Constant(rhs, shape=(1, 1)) 

1236 except Exception: 

1237 try: 

1238 rhs = Constant(rhs, shape=(lhs.shape[0], 1)) 

1239 except Exception: 

1240 rhs = Constant(rhs, shape=rhsShape) 

1241 else: 

1242 rhs = rhs.diag 

1243 else: 

1244 rhs = rhs.dupdiag(lhs.shape[0]) 

1245 elif rMatMul or lMatMul: 

1246 if lhs.shape == (1, 1): 

1247 # Any shape works. 

1248 rhs = Constant(rhs) 

1249 else: 

1250 try: # Try loading as scalar factor first. 

1251 rhs = Constant(rhs, shape=(1, 1)) 

1252 except Exception: 

1253 try: # Try loading as a vector. 

1254 if lMatMul: 

1255 s = (1, lhs.shape[0]) 

1256 else: 

1257 s = (lhs.shape[1], 1) 

1258 

1259 rhs = Constant(rhs, shape=s) 

1260 except TypeError: # Try loading as a matrix. 

1261 rhs = Constant(rhs, shape=rhsShape) 

1262 else: 

1263 rhs = Constant(rhs, shape=rhsShape) 

1264 except Exception as error: 

1265 fail("Could not load right hand side as a constant of " 

1266 "matching shape.", error) 

1267 elif rhsIsExpOrSet: 

1268 if allowNone and lhs is None: 

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

1270 

1271 lhsShape = None if rhsIsSet else make_lhs_shape(rhs.shape) 

1272 

1273 if diagBroadcast and rhs.shape[0] != rhs.shape[1]: 

1274 fail("Given that the left hand side operand is not a PICOS " 

1275 "expression, the right hand side operand must be square" 

1276 " for this operation.") 

1277 

1278 try: 

1279 if diagBroadcast: 

1280 try: 

1281 lhs = Constant(lhs, shape=(1, 1)) 

1282 except Exception: 

1283 try: 

1284 lhs = Constant(lhs, shape=(rhs.shape[0], 1)) 

1285 except Exception: 

1286 lhs = Constant(lhs, shape=lhsShape) 

1287 else: 

1288 lhs = lhs.diag 

1289 else: 

1290 lhs = lhs.dupdiag(rhs.shape[0]) 

1291 if rMatMul or lMatMul: 

1292 if rhs.shape == (1, 1): 

1293 # Any shape works. 

1294 lhs = Constant(lhs) 

1295 else: 

1296 try: # Try loading as scalar factor first. 

1297 lhs = Constant(lhs, shape=(1, 1)) 

1298 except TypeError: 

1299 try: # Try loading as a vector. 

1300 if rMatMul: 

1301 s = (1, rhs.shape[0]) 

1302 else: 

1303 s = (rhs.shape[1], 1) 

1304 

1305 lhs = Constant(lhs, shape=s) 

1306 except TypeError: # Try loading as a matrix. 

1307 lhs = Constant(lhs, shape=lhsShape) 

1308 else: 

1309 lhs = Constant(lhs, shape=lhsShape) 

1310 except Exception as error: 

1311 fail("Could not load left hand side as a constant of " 

1312 "matching shape.", error) 

1313 else: 

1314 assert False, "convert_operands is supposed to decorate " \ 

1315 "expression methods, but neither operand is a PICOS " \ 

1316 "expression or set." 

1317 

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

1319 return wrapper 

1320 return decorator 

1321 

1322 

1323def cvxopt_equals(A, B, absTol=None, relTol=None): 

1324 """Whether two CVXOPT (sparse) matrices are numerically equal or close. 

1325 

1326 For every common entry of ``A`` and ``B``, it is sufficient that one of the 

1327 two tolerances, ``absTol`` or ``relTol``, is satisfied. 

1328 

1329 :param float absTol: 

1330 Maximum allowed entrywise absolute difference. 

1331 

1332 :param float relTol: 

1333 Maximum allowed entrywise quotient of absolute difference at the entry 

1334 divided by the largest absolute value of any entry in both matrices. 

1335 """ 

1336 if A.size != B.size: 

1337 return False 

1338 

1339 Z = A - B 

1340 

1341 if not Z: 

1342 return True 

1343 

1344 if not absTol and not relTol: 

1345 return False 

1346 

1347 M = max(abs(Z)) 

1348 

1349 if relTol: 

1350 N = max(max(abs(A)), max(abs(B))) 

1351 

1352 if absTol and relTol: 

1353 if M > absTol and M / N > relTol: 

1354 return False 

1355 elif absTol: 

1356 if M > absTol: 

1357 return False 

1358 else: 

1359 if M / N > relTol: 

1360 return False 

1361 

1362 return True 

1363 

1364 

1365def cvxopt_maxdiff(A, B): 

1366 """Return the largest absolute difference of two (sparse) CVXOPT matrices. 

1367 

1368 :raises TypeError: If the matrices are not of the same shape. 

1369 """ 

1370 if A.size != B.size: 

1371 raise TypeError("The matrices do not have the same shape.") 

1372 

1373 # Work around "ValueError: max() arg is an empty sequence" for sparse zero. 

1374 if not A and not B: 

1375 return 0.0 

1376 

1377 return max(abs(A - B)) 

1378 

1379 

1380def cvxopt_hcat(matrices): 

1381 """Concatenate the given CVXOPT (sparse) matrices horizontally. 

1382 

1383 The resulting matrix is sparse if any input matrix is sparse. 

1384 

1385 :param list matrices: A list of CVXOPT (sparse) matrices. 

1386 """ 

1387 if not isinstance(matrices, list): 

1388 matrices = list(matrices) 

1389 

1390 sparse = any(isinstance(M, cvxopt.spmatrix) for M in matrices) 

1391 

1392 matrices = [[matrix] for matrix in matrices] 

1393 

1394 if sparse: 

1395 return cvxopt.sparse(matrices) 

1396 else: 

1397 return cvxopt.matrix(matrices) 

1398 

1399 

1400def cvxopt_vcat(matrices): 

1401 """Concatenate the given CVXOPT (sparse) matrices vertically. 

1402 

1403 The resulting matrix is sparse if any input matrix is sparse. 

1404 

1405 :param list matrices: A list of CVXOPT (sparse) matrices. 

1406 """ 

1407 if not isinstance(matrices, list): 

1408 matrices = list(matrices) 

1409 

1410 if any(isinstance(M, cvxopt.spmatrix) for M in matrices): 

1411 return cvxopt.sparse(matrices) 

1412 else: 

1413 return cvxopt.matrix(matrices) 

1414 

1415 

1416def cvxopt_hpsd(matrix): 

1417 """Whether the given CVXOPT matrix is hermitian positive semidefinite. 

1418 

1419 Uses :data:`~picos.settings.RELATIVE_HERMITIANNESS_TOLERANCE` and 

1420 :data:`~picos.settings.RELATIVE_SEMIDEFINITENESS_TOLERANCE`. 

1421 

1422 See also :func:`cvxopt_hpd`. 

1423 

1424 .. warning:: 

1425 

1426 The semidefiniteness tolerance allows negative, near-zero eigenvalues. 

1427 """ 

1428 if not cvxopt_equals(matrix, matrix.H, 

1429 relTol=settings.RELATIVE_HERMITIANNESS_TOLERANCE): 

1430 return False 

1431 

1432 eigenvalues = numpy.linalg.eigvalsh(cvx2np(matrix)) 

1433 

1434 minimum = -(max(list(eigenvalues) + [1.0]) 

1435 * settings.RELATIVE_SEMIDEFINITENESS_TOLERANCE) 

1436 

1437 return all(ev >= minimum for ev in eigenvalues) 

1438 

1439 

1440def cvxopt_hpd(matrix): 

1441 """Whether the given CVXOPT matrix is hermitian positive definite. 

1442 

1443 Uses :data:`~picos.settings.RELATIVE_HERMITIANNESS_TOLERANCE`. 

1444 

1445 See also :func:`cvxopt_hpsd`. 

1446 """ 

1447 if not cvxopt_equals(matrix, matrix.H, 

1448 relTol=settings.RELATIVE_HERMITIANNESS_TOLERANCE): 

1449 return False 

1450 

1451 eigenvalues = numpy.linalg.eigvalsh(cvx2np(matrix)) 

1452 

1453 return all(ev > 0 for ev in eigenvalues) 

1454 

1455 

1456def cvxopt_inverse(matrix): 

1457 """Return the inverse of the given CVXOPT matrix. 

1458 

1459 :raises ValueError: 

1460 If the matrix is not invertible. 

1461 """ 

1462 matrix_np = cvx2np(matrix) 

1463 

1464 try: 

1465 inverse_np = numpy.linalg.inv(matrix_np) 

1466 except numpy.linalg.LinAlgError as error: 

1467 raise ValueError("Failed to invert a {} CVXOPT matrix using NumPy." 

1468 .format(glyphs.shape(matrix.size))) from error 

1469 

1470 inverse, _ = load_data(inverse_np, matrix.size) 

1471 

1472 return inverse 

1473 

1474 

1475def cvxopt_principal_root(matrix): 

1476 """Return the principal square root of a symmetric positive semidef. matrix. 

1477 

1478 Given a real symmetric positive (semi)definite CVXOPT input matrix, returns 

1479 its unique positive (semi)definite matrix square root. 

1480 

1481 .. warning:: 

1482 

1483 Does not validate that the input matrix is symmetric positive 

1484 semidefinite and will still return a (useless) matrix if it is not. 

1485 """ 

1486 matrix_np = cvx2np(matrix) 

1487 U, s, _ = numpy.linalg.svd(matrix_np, hermitian=True) 

1488 root_np = numpy.dot(U*(s**0.5), U.T) 

1489 root, _ = load_data(root_np, matrix.size) 

1490 

1491 return root 

1492 

1493 

1494@lru_cache() 

1495def cvxopt_K(m, n, typecode="d"): 

1496 """The commutation matrix :math:`K_{(m,n)}` as a CVXOPT sparse matrix.""" 

1497 d = m*n 

1498 V = [1]*d 

1499 I = range(d) 

1500 J = [(k % n)*m + k // n for k in I] 

1501 return cvxopt.spmatrix(V, I, J, (d, d), typecode) 

1502 

1503 

1504def sparse_quadruple(A, reshape=None, preT=False, postT=False): 

1505 """Return a sparse representation of the given CVXOPT (sparse) matrix. 

1506 

1507 :param reshape: If set, then :math:`A` is reshaped on the fly. 

1508 :param bool preT: Transpose :math:`A` before reshaping. 

1509 :param bool postT: Transpose :math:`A` after reshaping. 

1510 

1511 :returns: A quadruple of values, row indices, column indices, and shape. 

1512 """ 

1513 if not isinstance(A, (cvxopt.spmatrix, cvxopt.matrix)): 

1514 raise TypeError("Input must be a CVXOPT (sparse) matrix.") 

1515 

1516 m, n = A.size 

1517 

1518 if reshape: 

1519 reshape = load_shape(reshape) 

1520 

1521 if A.size[0] * A.size[1] != reshape[0] * reshape[1]: 

1522 raise TypeError("Cannot reshape from {} to {}.".format( 

1523 glyphs.shape(A.size), glyphs.shape(reshape))) 

1524 

1525 if not A: 

1526 V, I, J = [], [], [] 

1527 p, q = reshape if reshape else (n, m if preT else m, n) 

1528 else: 

1529 if isinstance(A, cvxopt.matrix): 

1530 A = cvxopt.sparse(A) 

1531 

1532 V, I, J = A.V, A.I, A.J 

1533 

1534 if preT: 

1535 I, J, m, n = J, I, n, m 

1536 

1537 if reshape: 

1538 p, q = reshape 

1539 I, J = zip(*[((i + j*m) % p, (i + j*m) // p) for i, j in zip(I, J)]) 

1540 else: 

1541 p, q = m, n 

1542 

1543 if postT: 

1544 I, J, p, q = J, I, q, p 

1545 

1546 return V, I, J, (p, q) 

1547 

1548 

1549def left_kronecker_I(A, k, reshape=None, preT=False, postT=False): 

1550 r"""Return :math:`I_k \otimes A` for a CVXOPT (sparse) matrix :math:`A`. 

1551 

1552 In other words, if :math:`A` is a :math:`m \times n` CVXOPT (sparse) 

1553 matrix, returns a :math:`km \times kn` CVXOPT sparse block matrix with 

1554 all blocks of size :math:`m \times n`, the diagonal blocks (horizontal 

1555 block index equal to vertical block index) equal to :math:`A`, and all 

1556 other blocks zero. 

1557 

1558 :param reshape: If set, then :math:`A` is reshaped on the fly. 

1559 :param bool preT: Transpose :math:`A` before reshaping. 

1560 :param bool postT: Transpose :math:`A` after reshaping. 

1561 

1562 :returns: 

1563 If :math:`A` is dense and :math:`k = 1`, a 

1564 :func:`CVXOPT dense matrix <cvxopt:cvxopt.matrix>`, otherwise a 

1565 :func:`CVXOPT sparse matrix <cvxopt:cvxopt.spmatrix>`. 

1566 """ 

1567 A = A.T if preT else A 

1568 

1569 if reshape: 

1570 A = A if preT else A[:] # Copy if not already fresh. 

1571 A.size = reshape 

1572 

1573 A = A.T if postT else A 

1574 

1575 if k == 1: 

1576 if not preT and not reshape and not postT: 

1577 return A + 0 # Always return a copy. 

1578 else: 

1579 return A 

1580 

1581 if A.size[0] == A.size[1]: 

1582 A = cvxopt.spdiag([A for _ in range(k)]) 

1583 else: 

1584 Z = cvxopt.spmatrix([], [], [], A.size) 

1585 A = cvxopt.sparse([[Z]*i + [A] + [Z]*(k-i-1) for i in range(k)]) 

1586 

1587 return A.T if postT else A 

1588 

1589 

1590def right_kronecker_I(A, k, reshape=None, preT=False, postT=False): 

1591 r"""Return :math:`A \otimes I_k` for a CVXOPT (sparse) matrix :math:`A`. 

1592 

1593 :param reshape: If set, then :math:`A` is reshaped on the fly. 

1594 :param bool preT: Transpose :math:`A` before reshaping. 

1595 :param bool postT: Transpose :math:`A` after reshaping. 

1596 

1597 :returns: 

1598 If :math:`A` is dense and :math:`k = 1`, a 

1599 :func:`CVXOPT dense matrix <cvxopt:cvxopt.matrix>`, otherwise a 

1600 :func:`CVXOPT sparse matrix <cvxopt:cvxopt.spmatrix>`. 

1601 """ 

1602 if isinstance(A, cvxopt.matrix): 

1603 # Dense case: Use NumPy. 

1604 A = numpy.array(A) 

1605 

1606 A = A.T if preT else A 

1607 A = A.reshape(reshape, order="F") if reshape else A 

1608 A = A.T if postT else A 

1609 

1610 A = numpy.kron(A, numpy.eye(k)) 

1611 

1612 return load_data(A, sparse=(k > 1))[0] 

1613 else: 

1614 # Sparse case: Python implementation. 

1615 # This is slower than the NumPy approach in general but can handle 

1616 # arbitrarily large matrices given that they are sufficiently sparse. 

1617 V, I, J, shape = sparse_quadruple(A, reshape, preT, postT) 

1618 m, n = shape 

1619 

1620 if V: 

1621 V, I, J = zip(*[(v, k*i + l, k*j + l) 

1622 for v, i, j in zip(V, I, J) for l in range(k)]) 

1623 

1624 return cvxopt.spmatrix(V, I, J, (k*m, k*n)) 

1625 

1626 

1627def cvx2np(A, reshape=None): 

1628 """Convert a CVXOPT (sparse) matrix to a NumPy two-dimensional array. 

1629 

1630 :param A: The CVXOPT :func:`dense <cvxopt:cvxopt.matrix>` or 

1631 :func:`sparse <cvxopt:cvxopt.spmatrix>` matrix to convert. 

1632 :param bool reshape: Optional new shape for the converted matrix. 

1633 

1634 :returns: Converted :class:`NumPy array <numpy.ndarray>`. 

1635 """ 

1636 assert isinstance(A, (cvxopt.matrix, cvxopt.spmatrix)) 

1637 

1638 if isinstance(A, cvxopt.spmatrix): 

1639 A = cvxopt.matrix(A) 

1640 

1641 if reshape: 

1642 shape = load_shape(reshape) 

1643 return numpy.reshape(A, shape, "F") 

1644 else: 

1645 return numpy.array(A) 

1646 

1647 

1648def cvx2csc(A): 

1649 """Convert a CVXOPT matrix to a SciPy sparse matrix in CSC format.""" 

1650 import scipy.sparse 

1651 

1652 assert isinstance(A, (cvxopt.matrix, cvxopt.spmatrix)) 

1653 

1654 if isinstance(A, cvxopt.spmatrix): 

1655 csc = tuple(tuple(x) for x in reversed(A.CCS)) 

1656 return scipy.sparse.csc_matrix(csc, shape=A.size) 

1657 else: 

1658 return scipy.sparse.csc_matrix(A) 

1659 

1660 

1661def cvx2csr(A): 

1662 """Convert a CVXOPT matrix to a SciPy sparse matrix in CSR format.""" 

1663 import scipy.sparse 

1664 

1665 assert isinstance(A, (cvxopt.matrix, cvxopt.spmatrix)) 

1666 

1667 if isinstance(A, cvxopt.spmatrix): 

1668 csc = tuple(tuple(x) for x in reversed(A.CCS)) 

1669 csc = scipy.sparse.csc_matrix(csc, shape=A.size) 

1670 return csc.tocsr() 

1671 else: 

1672 return scipy.sparse.csr_matrix(A) 

1673 

1674 

1675def sp2cvx(A): 

1676 """Convert a SciPy sparse matrix to a CVXOPT sparse matrix.""" 

1677 import scipy.sparse 

1678 

1679 assert isinstance(A, scipy.sparse.spmatrix) 

1680 

1681 A = A.tocoo() 

1682 

1683 V = A.data 

1684 I = A.row 

1685 J = A.col 

1686 

1687 if issubclass(A.dtype.type, numpy.complex): 

1688 tc = "z" 

1689 else: 

1690 tc = "d" 

1691 

1692 return cvxopt.spmatrix(V, I, J, A.shape, tc) 

1693 

1694 

1695def make_fraction(p, denominator_limit): 

1696 """Convert a float :math:`p` to a limited precision fraction. 

1697 

1698 :param float p: The float to convert, may be positive or negative infinity. 

1699 :param int denominator_limit: The largest allowed denominator. 

1700 

1701 :returns tuple: A quadruple ``(num, den, pNew, pStr)`` with ``pNew`` the 

1702 limited precision version of :math:`p`, ``pStr`` a string representation 

1703 of the fraction, and ``num`` and ``den`` the numerator and the 

1704 denominator of the fraction, respectively. 

1705 """ 

1706 # LEGACY: Old tools.tracepow allowed this. 

1707 if p in ("inf", "-inf"): 

1708 p = float(p) 

1709 

1710 if p in (float("inf"), float("-inf")): 

1711 return p, 1, p, glyphs.neg(glyphs.infty) if p < 0 else glyphs.infty 

1712 

1713 frac = Fraction(p).limit_denominator(denominator_limit) 

1714 num = frac.numerator 

1715 den = frac.denominator 

1716 pNew = float(num) / float(den) 

1717 

1718 if den == 1: 

1719 pStr = glyphs.scalar(num) 

1720 else: 

1721 pStr = glyphs.clever_div(glyphs.scalar(num), glyphs.scalar(den)) 

1722 

1723 return num, den, pNew, pStr 

1724 

1725 

1726def value(obj, sparse=None, numpy=False): 

1727 """Convert (nested) PICOS objects to their current value. 

1728 

1729 :param obj: 

1730 Either a single (PICOS) object that has a ``value`` attribute, such as a 

1731 :class:`mutable <picos.expressions.mutable.Mutable>`, 

1732 :class:`expression <picos.expressions.expression.Expression>` or 

1733 :class:`~picos.modeling.Problem`, or a (nested) :class:`list`, 

1734 :class:`tuple` or :class:`dict` thereof. 

1735 

1736 :param sparse: 

1737 If :obj:`None`, retrieved multidimensional values can be returned as 

1738 either CVXOPT :func:`sparse <cvxopt:cvxopt.spmatrix>` or 

1739 :func:`dense <cvxopt:cvxopt.matrix>` matrices, whichever PICOS stores 

1740 internally. If :obj:`True` or :obj:`False`, multidimensional values are 

1741 always returned as sparse or dense types, respectively. 

1742 

1743 :param bool numpy: 

1744 If :obj:`True`, retrieved multidimensional values are returned as a 

1745 NumPy :class:`~numpy:numpy.ndarray` instead of a CVXOPT type. May not be 

1746 set in combination with ``sparse=True``. 

1747 

1748 :returns: 

1749 An object of the same (nested) structure as ``obj``, with every 

1750 occurence of any object with a ``value`` attribute replaced by that 

1751 attribute's current numeric value. In the case of dictionaries, only the 

1752 dictionary values will be converted. 

1753 

1754 :raises TypeError: 

1755 If some object with a ``value`` attribute has a value that cannot be 

1756 converted to a matrix by :func:`~.load_data`. This can only happen if 

1757 the object in question is not a PICOS object. 

1758 

1759 :Example: 

1760 

1761 >>> from picos import RealVariable, value 

1762 >>> from pprint import pprint 

1763 >>> x = {key: RealVariable(key) for key in ("foo", "bar")} 

1764 >>> x["foo"].value = 2 

1765 >>> x["bar"].value = 3 

1766 >>> pprint(value(x)) 

1767 {'bar': 3.0, 'foo': 2.0} 

1768 """ 

1769 if sparse and numpy: 

1770 raise ValueError("NumPy does not support sparse matrices.") 

1771 

1772 if isinstance(obj, tuple): 

1773 return tuple(value(inner, sparse, numpy) for inner in obj) 

1774 elif isinstance(obj, list): 

1775 return [value(inner, sparse, numpy) for inner in obj] 

1776 elif isinstance(obj, dict): 

1777 return {k: value(v, sparse, numpy) for k, v in obj.items()} 

1778 else: 

1779 if hasattr(obj, "value"): 

1780 val = obj.value 

1781 

1782 if isinstance(val, (int, float, complex)): 

1783 return val 

1784 

1785 # PICOS objects always return their value as a CVXOPT matrix type, 

1786 # but this function may be used on other objects with a value 

1787 # attribute. Try to convert their value to a CVXOPT matrix first. 

1788 if not isinstance(val, (cvxopt.matrix, cvxopt.spmatrix)): 

1789 if numpy: 

1790 load_sparse = False 

1791 elif sparse: 

1792 load_sparse = True 

1793 else: 

1794 load_sparse = None 

1795 

1796 val = load_data(val, sparse=load_sparse)[0] 

1797 elif isinstance(val, cvxopt.spmatrix) \ 

1798 and (sparse is False or numpy): 

1799 val = cvxopt.dense(val) 

1800 

1801 if numpy: 

1802 import numpy as the_numpy 

1803 assert isinstance(val, cvxopt.matrix) 

1804 val = the_numpy.array(val) 

1805 

1806 return val 

1807 else: 

1808 return obj 

1809 

1810 

1811# -------------------------------------- 

1812__all__ = api_end(_API_START, globals())