Hide keyboard shortcuts

Hot-keys on this page

r m x p   toggle line displays

j k   next/prev highlighted chunk

0   (zero) top of page

1   (one) first highlighted chunk

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

2# Copyright (C) 2019-2020 Maximilian Stahlberg 

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

4# 

5# This file is part of PICOS. 

6# 

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

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

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

10# version. 

11# 

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

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

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

15# 

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

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

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

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 load_data(value, shape=None, typecode=None, sparse=None, 

129 alwaysCopy=True, legacy=False): 

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

131 

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

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

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

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

136 worrying about any conversion. 

137 

138 :Supported data types: 

139 

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: 

150 

151 .. code-block:: python 

152 

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

154 [4, 5, 6]] 

155 

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: 

159 

160 .. code-block:: python 

161 

162 A = '''1 2 3 

163 4 5 6''' 

164 

165 - An algebraic string description: 

166 

167 .. list-table:: 

168 :widths: 1 99 

169 

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

187 

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. 

194 

195 :Broadcasting and reshaping rules: 

196 

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

219 

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 

227 

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. 

233 

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 

239 

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. 

244 

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. 

249 

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. 

253 

254 :Example: 

255 

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

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

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

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

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

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

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

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

264 >>> load_data(A) 

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

266 >>> # Data as string: 

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

268 >>> print(string) 

269 [7×2:e_7,2] 

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

271 [ 0 0 ] 

272 [ 0 0 ] 

273 [ 0 0 ] 

274 [ 0 0 ] 

275 [ 0 0 ] 

276 [ 0 0 ] 

277 [ 0 1.00e+00] 

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

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

280 """ 

281 from .exp_affine import ComplexAffineExpression 

282 from .expression import Expression 

283 

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 

290 

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 

306 

307 raise TypeError( 

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

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

310 

311 def load_dense(value, shape, typecode): 

312 """Create a CVXOPT dense matrix.""" 

313 try: 

314 if typecode: 

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

316 else: 

317 return cvxopt.matrix(value, shape) 

318 except TypeError as error: 

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

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

321 if typecode in "id": 

322 try: 

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

324 except Exception: 

325 pass 

326 else: 

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

328 if typecode == "d": 

329 return complexValue.real() 

330 else: 

331 realData = list(complexValue.real()) 

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

333 if intData == realData: 

334 try: 

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

336 except Exception: 

337 pass 

338 

339 raise TypeError( 

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

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

342 

343 def simple_vector_as_row(shape): 

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

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

346 

347 def broadcast_error(): 

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

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

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

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

352 

353 def scalar(x, typecode=None): 

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

355 x = complex(x) 

356 

357 if typecode == "z": 

358 return x # We are content with complex. 

359 else: 

360 if not x.imag: 

361 x = x.real 

362 

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. 

372 

373 raise TypeError( 

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

375 

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

377 # where None means any dimensionality. 

378 shape = load_shape(shape, wildcards=True) 

379 

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

384 

385 # Validate the sparsity argument. 

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

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

388 

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 

396 

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

398 inShape = None 

399 

400 # Allow conversions to specify their own string description. 

401 string = None 

402 

403 # Convert from range to list. 

404 if isinstance(value, range): 

405 value = list(value) 

406 

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

408 if isinstance(value, Expression): 

409 value = value.refined 

410 

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

424 

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

430 

431 inShape = (1, 1) 

432 outShape = blend_shapes(shape, inShape) 

433 

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) 

437 

438 if not value and sparse is not False: 

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

440 else: 

441 value = load_dense(value, outShape, typecode) 

442 

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. 

452 

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) 

457 

458 # Refine the output shape. 

459 inShape = value.size 

460 outShape = blend_shapes(shape, inShape) 

461 

462 # Define shorthands. 

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

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

465 

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 

471 

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

473 reshaped = True 

474 if inShape == outShape: 

475 reshaped = False 

476 elif inShape == (1, 1): 

477 # Broadcast a scalar. 

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

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

480 # Copy columns horizontally. 

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

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

483 # Copy rows vertically. 

484 value = load_dense( 

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

486 elif (inLength, 1) == outShape: 

487 # Vectorize in column-major order. 

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

489 # Quick method. 

490 value = value[:] 

491 else: 

492 # With typecode change. 

493 value = load_dense(value, outShape, typecode) 

494 elif (1, inLength) == outShape: 

495 # Vectorize in row-major order. 

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

497 # Quick method. 

498 value = value.T[:].T 

499 else: 

500 # With typecode change. 

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

502 elif (outLength, 1) == inShape: 

503 # Devectorize in column-major order. 

504 value = load_dense(value, outShape, typecode) 

505 elif (1, outLength) == inShape: 

506 # Devectorize in row-major order. 

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

508 else: 

509 broadcast_error() 

510 if reshaped: 

511 alwaysCopy = False 

512 

513 # The data now has the desired shape. 

514 assert value.size == outShape 

515 

516 # Ensure the proper typecode. 

517 if typecode and value.typecode != typecode: 

518 value = load_dense(value, outShape, typecode) 

519 alwaysCopy = False 

520 

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. 

526 

527 # Refine the output shape. 

528 inShape = value.size 

529 outShape = blend_shapes(shape, inShape) 

530 

531 # Define shorthands. 

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

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

534 

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 

540 

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

542 reshaped = True 

543 if inShape == outShape: 

544 reshaped = False 

545 elif inShape == (1, 1): 

546 # Broadcast a scalar. 

547 if value[0] != 0: 

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

549 else: 

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

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

552 # Copy columns horizontally. 

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

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

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

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

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

558 # Copy rows vertically. 

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

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

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

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

563 elif (inLength, 1) == outShape: 

564 # Vectorize in column-major order. 

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

566 # Quick method. 

567 value = value[:] 

568 else: 

569 # With typecode change. 

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

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

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

573 J = [0]*nnz 

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

575 elif (1, inLength) == outShape: 

576 # Vectorize in row-major order. 

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

578 # Quick method. 

579 value = value.T[:].T 

580 else: 

581 # With typecode change. 

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

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

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

585 I = [0]*nnz 

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

587 elif (outLength, 1) == inShape: 

588 # Devectorize in column-major order. 

589 # TODO: Logic for also changing the typecode. 

590 value = value[:] 

591 value.size = outShape 

592 elif (1, outLength) == inShape: 

593 # Devectorize in row-major order. 

594 # TODO: Logic for also changing the typecode. 

595 value = value[:] 

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

597 value = value.T 

598 else: 

599 broadcast_error() 

600 if reshaped: 

601 alwaysCopy = False 

602 

603 # The data now has the desired shape. 

604 assert value.size == outShape 

605 

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 

613 

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

629 

630 # Refine the output shape. 

631 inShape = load_shape(value.shape) 

632 outShape = blend_shapes(shape, inShape) 

633 

634 # Define shorthands. 

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

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

637 

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) 

644 

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 

650 

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

652 if inShape == outShape: 

653 pass 

654 elif inShape == (1, 1): 

655 # Broadcast a scalar. 

656 value = numpy.full(outShape, value) 

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

658 # Copy columns horizontally. 

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

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

661 # Copy rows vertically. 

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

663 elif (inLength, 1) == outShape: 

664 # Vectorize in column-major order. 

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

666 elif (1, inLength) == outShape: 

667 # Vectorize in row-major order. 

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

669 elif (outLength, 1) == inShape: 

670 # Devectorize in column-major order. 

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

672 elif (1, outLength) == inShape: 

673 # Devectorize in row-major order. 

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

675 else: 

676 broadcast_error() 

677 

678 # The data now has the desired shape. 

679 assert value.shape == outShape 

680 

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 

684 

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

697 

698 # Retrieve a copy of the numeric value. 

699 value = value.value 

700 

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] 

707 

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) 

718 

719 match = _LOAD_DATA_REGEX.match(value) 

720 

721 if not match: 

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

723 .format(value)) 

724 

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 

730 

731 # Convert the factor. 

732 factor = scalar(factor) if factor else 1 

733 

734 # Determine whether the matrix will be square. 

735 square = base == "I" 

736 

737 # Convert the shape. 

738 if inShape: 

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

740 if len(inShape) == 1: 

741 inShape *= 2 

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

743 

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

748 

749 outShape = inShape 

750 elif None in shape: 

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

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

753 assert None not in outShape 

754 else: 

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

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

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

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

759 else: 

760 outShape = shape 

761 

762 # Create the base matrix. 

763 if base.startswith("e_"): 

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

765 if len(position) == 1: 

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

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

768 else: 

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

770 

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

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

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

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

775 

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

777 

778 if outShape[1] == 1: 

779 assert j == 0 

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

781 elif outShape[0] == 1: 

782 assert i == 0 

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

784 else: 

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

786 elif base.startswith("|"): 

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

788 

789 # Pull the factor inside the matrix. 

790 element *= factor 

791 factor = 1 

792 

793 value = load_dense(element, outShape, typecode) 

794 string = glyphs.scalar(element) 

795 elif base == "I": 

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

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

798 

799 n = outShape[0] 

800 

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) 

805 

806 # Apply a coefficient. 

807 if factor != 1: 

808 value = value * factor 

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

810 

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

821 

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

830 

831 rows = len(value) 

832 cols = None 

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

834 

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) 

839 

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

845 

846 # Make sure row length is consistent. 

847 if cols is None: 

848 cols = len(row) 

849 

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

856 

857 if isinstance(row, tuple): 

858 value[rowNum] = list(row) 

859 

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

861 elif rows == 1: 

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

863 return load_data( 

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

865 else: 

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

867 value = load_dense(value, outShape, typecode) 

868 

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

874 

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) 

882 

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 

889 

890 value = cvxopt.spdiag(value) 

891 

892 if sparse is False: 

893 value = cvxopt.dense(value) 

894 

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

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

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

898 

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) 

908 

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

910 if string is None: 

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

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

913 else: 

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

915 

916 return value, string 

917 

918 

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

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

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

922 alwaysCopy=alwaysCopy) 

923 

924 

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

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

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

928 alwaysCopy=alwaysCopy) 

929 

930 

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

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

933 

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

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

936 made to load them as constant expressions. 

937 

938 :Decorator guarantee: 

939 

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

941 exectued. 

942 

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 

953 

954 def name(): 

955 if hasattr(func, "__qualname__"): 

956 return func.__qualname__ 

957 else: 

958 return func.__name__ 

959 

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

961 

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 

978 

979 newargs[key] = arg 

980 

981 return func(**newargs) 

982 return wrapper 

983 return decorator 

984 

985 

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. 

990 

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. 

995 

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. 

1001 

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. 

1007 

1008 :Decorator guarantee: 

1009 

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. 

1014 

1015 :Broadcasting rules for affine expressions: 

1016 

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. 

1021 

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. 

1034 

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

1036 operands being raw data or an affine expression. 

1037 

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 

1047 

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

1054 

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 

1065 

1066 def make_lhs_shape(rhsShape): 

1067 if sameShape: 

1068 return rhsShape 

1069 elif horiCat: 

1070 return (rhsShape[0], None) 

1071 elif vertCat: 

1072 return (None, rhsShape[1]) 

1073 elif rMatMul: 

1074 return (None, rhsShape[0]) 

1075 elif lMatMul: 

1076 return (rhsShape[1], None) 

1077 elif scalarLHS: 

1078 return (1, 1) 

1079 else: 

1080 return None 

1081 

1082 def make_rhs_shape(lhsShape): 

1083 if sameShape: 

1084 return lhsShape 

1085 elif horiCat: 

1086 return (lhsShape[0], None) 

1087 elif vertCat: 

1088 return (None, lhsShape[1]) 

1089 elif rMatMul: 

1090 return (lhsShape[1], None) 

1091 elif lMatMul: 

1092 return (None, lhsShape[0]) 

1093 elif scalarRHS: 

1094 return (1, 1) 

1095 else: 

1096 return None 

1097 

1098 from .exp_affine import Constant 

1099 from .exp_biaffine import BiaffineExpression 

1100 from .expression import Expression 

1101 from .set import Set 

1102 

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

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

1105 

1106 lhsIsSet = lhsIsExpOrSet and isinstance(lhs, Set) 

1107 rhsIsSet = rhsIsExpOrSet and isinstance(rhs, Set) 

1108 

1109 if lhsIsSet: 

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

1111 "not pose any shape requirements." 

1112 

1113 if lhsIsExpOrSet and rhsIsExpOrSet: 

1114 lhsIsAffine = isinstance(lhs, BiaffineExpression) 

1115 rhsIsAffine = isinstance(rhs, BiaffineExpression) 

1116 

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) 

1121 

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

1123 # already expressions. 

1124 if not selected: 

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

1126 

1127 assert not lhsIsSet # Handled earlier. 

1128 

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

1130 if rhsIsSet: 

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

1132 

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

1134 

1135 # Check if already matching size. 

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

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

1138 

1139 lm, ln = lhs.shape 

1140 rm, rn = rhs.shape 

1141 

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) 

1148 

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

1150 

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) 

1160 

1161 # Matrix multiplication always accepts scalars. 

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

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

1164 

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

1166 if lhsIsAffine and lhsL == 1: 

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

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

1169 

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

1171 if rhsIsAffine and rhsL == 1: 

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

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

1174 

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) 

1182 

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

1184 

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) 

1203 

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

1205 

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

1225 

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

1227 return wrapper 

1228 return decorator 

1229 

1230 

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

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

1233 

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

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

1236 

1237 :param float absTol: 

1238 Maximum allowed entrywise absolute difference. 

1239 

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 

1246 

1247 Z = A - B 

1248 

1249 if not Z: 

1250 return True 

1251 

1252 if not absTol and not relTol: 

1253 return False 

1254 

1255 M = max(abs(Z)) 

1256 

1257 if relTol: 

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

1259 

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 

1269 

1270 return True 

1271 

1272 

1273def cvxopt_maxdiff(A, B): 

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

1275 

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

1280 

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

1282 if not A and not B: 

1283 return 0.0 

1284 

1285 return max(abs(A - B)) 

1286 

1287 

1288def cvxopt_hcat(matrices): 

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

1290 

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

1292 

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

1294 """ 

1295 if not isinstance(matrices, list): 

1296 matrices = list(matrices) 

1297 

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

1299 

1300 matrices = [[matrix] for matrix in matrices] 

1301 

1302 if sparse: 

1303 return cvxopt.sparse(matrices) 

1304 else: 

1305 return cvxopt.matrix(matrices) 

1306 

1307 

1308def cvxopt_vcat(matrices): 

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

1310 

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

1312 

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

1314 """ 

1315 if not isinstance(matrices, list): 

1316 matrices = list(matrices) 

1317 

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

1319 return cvxopt.sparse(matrices) 

1320 else: 

1321 return cvxopt.matrix(matrices) 

1322 

1323 

1324def cvxopt_hpsd(matrix): 

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

1326 

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

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

1329 

1330 See also :func:`cvxopt_hpd`. 

1331 

1332 .. warning:: 

1333 

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 

1339 

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

1341 

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

1343 * settings.RELATIVE_SEMIDEFINITENESS_TOLERANCE) 

1344 

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

1346 

1347 

1348def cvxopt_hpd(matrix): 

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

1350 

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

1352 

1353 See also :func:`cvxopt_hpsd`. 

1354 """ 

1355 if not cvxopt_equals(matrix, matrix.H, 

1356 relTol=settings.RELATIVE_HERMITIANNESS_TOLERANCE): 

1357 return False 

1358 

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

1360 

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

1362 

1363 

1364def cvxopt_inverse(matrix): 

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

1366 

1367 :raises ValueError: 

1368 If the matrix is not invertible. 

1369 """ 

1370 matrix_np = cvx2np(matrix) 

1371 

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 

1377 

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

1379 

1380 return inverse 

1381 

1382 

1383def cvxopt_principal_root(matrix): 

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

1385 

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

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

1388 

1389 .. warning:: 

1390 

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) 

1398 

1399 return root 

1400 

1401 

1402@lru_cache() 

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

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

1405 d = m*n 

1406 V = [1]*d 

1407 I = range(d) 

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

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

1410 

1411 

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

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

1414 

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. 

1418 

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

1423 

1424 m, n = A.size 

1425 

1426 if reshape: 

1427 reshape = load_shape(reshape) 

1428 

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

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

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

1432 

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 

1438 

1439 if preT: 

1440 I, J, m, n = J, I, n, m 

1441 

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 

1447 

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

1468 

1469 if postT: 

1470 I, J, p, q = J, I, q, p 

1471 

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

1473 

1474 

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

1477 

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. 

1483 

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. 

1487 

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 

1494 

1495 if reshape: 

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

1497 A.size = reshape 

1498 

1499 A = A.T if postT else A 

1500 

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 

1506 

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

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

1509 else: 

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

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

1512 

1513 return A.T if postT else A 

1514 

1515 

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

1518 

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. 

1522 

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) 

1531 

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 

1535 

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

1537 

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

1539 else: 

1540 # Sparse case: Python implementation. 

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

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

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

1544 m, n = shape 

1545 

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

1549 

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

1551 

1552 

1553def cvx2np(A, reshape=None): 

1554 """Convert a CVXOPT (sparse) matrix to a NumPy matrix. 

1555 

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. 

1559 

1560 :returns: Converted :class:`NumPy matrix <numpy.matrix>`. 

1561 """ 

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

1563 

1564 if isinstance(A, cvxopt.spmatrix): 

1565 A = cvxopt.matrix(A) 

1566 

1567 if reshape: 

1568 shape = load_shape(reshape) 

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

1570 else: 

1571 return numpy.matrix(A) 

1572 

1573 

1574def make_fraction(p, denominator_limit): 

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

1576 

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

1578 :param int denominator_limit: The largest allowed denominator. 

1579 

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) 

1588 

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

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

1591 

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

1593 num = frac.numerator 

1594 den = frac.denominator 

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

1596 

1597 if den == 1: 

1598 pStr = glyphs.scalar(num) 

1599 else: 

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

1601 

1602 return num, den, pNew, pStr 

1603 

1604 

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

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

1607 

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. 

1614 

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. 

1621 

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

1626 

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. 

1632 

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. 

1637 

1638 :Example: 

1639 

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

1650 

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 

1660 

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

1662 return val 

1663 

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

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

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

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

1668 if numpy: 

1669 load_sparse = False 

1670 elif sparse: 

1671 load_sparse = True 

1672 else: 

1673 load_sparse = None 

1674 

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

1676 elif isinstance(val, cvxopt.spmatrix) \ 

1677 and (sparse is False or numpy): 

1678 val = cvxopt.dense(val) 

1679 

1680 if numpy: 

1681 import numpy as the_numpy 

1682 assert isinstance(val, cvxopt.matrix) 

1683 val = the_numpy.array(val) 

1684 

1685 return val 

1686 else: 

1687 return obj 

1688 

1689 

1690# -------------------------------------- 

1691__all__ = api_end(_API_START, globals())