Coverage for /usr/lib/python3/dist-packages/mpmath/matrices/matrices.py: 15%

1from ..libmp.backend import xrange 

2import warnings 

3 

4# TODO: interpret list as vectors (for multiplication) 

5 

6rowsep = '\n' 

7colsep = ' ' 

8 

9class _matrix(object): 

10 """ 

11 Numerical matrix. 

12 

13 Specify the dimensions or the data as a nested list. 

14 Elements default to zero. 

15 Use a flat list to create a column vector easily. 

16 

17 The datatype of the context (mpf for mp, mpi for iv, and float for fp) is used to store the data. 

18 

19 Creating matrices 

20 ----------------- 

21 

22 Matrices in mpmath are implemented using dictionaries. Only non-zero values 

23 are stored, so it is cheap to represent sparse matrices. 

24 

25 The most basic way to create one is to use the ``matrix`` class directly. 

26 You can create an empty matrix specifying the dimensions: 

27 

28 >>> from mpmath import * 

29 >>> mp.dps = 15 

30 >>> matrix(2) 

31 matrix( 

32 [['0.0', '0.0'], 

33 ['0.0', '0.0']]) 

34 >>> matrix(2, 3) 

35 matrix( 

36 [['0.0', '0.0', '0.0'], 

37 ['0.0', '0.0', '0.0']]) 

38 

39 Calling ``matrix`` with one dimension will create a square matrix. 

40 

41 To access the dimensions of a matrix, use the ``rows`` or ``cols`` keyword: 

42 

43 >>> A = matrix(3, 2) 

44 >>> A 

45 matrix( 

46 [['0.0', '0.0'], 

47 ['0.0', '0.0'], 

48 ['0.0', '0.0']]) 

49 >>> A.rows 

50 3 

51 >>> A.cols 

52 2 

53 

54 You can also change the dimension of an existing matrix. This will set the 

55 new elements to 0. If the new dimension is smaller than before, the 

56 concerning elements are discarded: 

57 

58 >>> A.rows = 2 

59 >>> A 

60 matrix( 

61 [['0.0', '0.0'], 

62 ['0.0', '0.0']]) 

63 

64 Internally ``mpmathify`` is used every time an element is set. This 

65 is done using the syntax A[row,column], counting from 0: 

66 

67 >>> A = matrix(2) 

68 >>> A[1,1] = 1 + 1j 

69 >>> A 

70 matrix( 

71 [['0.0', '0.0'], 

72 ['0.0', mpc(real='1.0', imag='1.0')]]) 

73 

74 A more comfortable way to create a matrix lets you use nested lists: 

75 

76 >>> matrix([[1, 2], [3, 4]]) 

77 matrix( 

78 [['1.0', '2.0'], 

79 ['3.0', '4.0']]) 

80 

81 Convenient advanced functions are available for creating various standard 

82 matrices, see ``zeros``, ``ones``, ``diag``, ``eye``, ``randmatrix`` and 

83 ``hilbert``. 

84 

85 Vectors 

86 ....... 

87 

88 Vectors may also be represented by the ``matrix`` class (with rows = 1 or cols = 1). 

89 For vectors there are some things which make life easier. A column vector can 

90 be created using a flat list, a row vectors using an almost flat nested list:: 

91 

92 >>> matrix([1, 2, 3]) 

93 matrix( 

94 [['1.0'], 

95 ['2.0'], 

96 ['3.0']]) 

97 >>> matrix([[1, 2, 3]]) 

98 matrix( 

99 [['1.0', '2.0', '3.0']]) 

100 

101 Optionally vectors can be accessed like lists, using only a single index:: 

102 

103 >>> x = matrix([1, 2, 3]) 

104 >>> x[1] 

105 mpf('2.0') 

106 >>> x[1,0] 

107 mpf('2.0') 

108 

109 Other 

110 ..... 

111 

112 Like you probably expected, matrices can be printed:: 

113 

114 >>> print randmatrix(3) # doctest:+SKIP 

115 [ 0.782963853573023 0.802057689719883 0.427895717335467] 

116 [0.0541876859348597 0.708243266653103 0.615134039977379] 

117 [ 0.856151514955773 0.544759264818486 0.686210904770947] 

118 

119 Use ``nstr`` or ``nprint`` to specify the number of digits to print:: 

120 

121 >>> nprint(randmatrix(5), 3) # doctest:+SKIP 

122 [2.07e-1 1.66e-1 5.06e-1 1.89e-1 8.29e-1] 

123 [6.62e-1 6.55e-1 4.47e-1 4.82e-1 2.06e-2] 

124 [4.33e-1 7.75e-1 6.93e-2 2.86e-1 5.71e-1] 

125 [1.01e-1 2.53e-1 6.13e-1 3.32e-1 2.59e-1] 

126 [1.56e-1 7.27e-2 6.05e-1 6.67e-2 2.79e-1] 

127 

128 As matrices are mutable, you will need to copy them sometimes:: 

129 

130 >>> A = matrix(2) 

131 >>> A 

132 matrix( 

133 [['0.0', '0.0'], 

134 ['0.0', '0.0']]) 

135 >>> B = A.copy() 

136 >>> B[0,0] = 1 

137 >>> B 

138 matrix( 

139 [['1.0', '0.0'], 

140 ['0.0', '0.0']]) 

141 >>> A 

142 matrix( 

143 [['0.0', '0.0'], 

144 ['0.0', '0.0']]) 

145 

146 Finally, it is possible to convert a matrix to a nested list. This is very useful, 

147 as most Python libraries involving matrices or arrays (namely NumPy or SymPy) 

148 support this format:: 

149 

150 >>> B.tolist() 

151 [[mpf('1.0'), mpf('0.0')], [mpf('0.0'), mpf('0.0')]] 

152 

153 

154 Matrix operations 

155 ----------------- 

156 

157 You can add and subtract matrices of compatible dimensions:: 

158 

159 >>> A = matrix([[1, 2], [3, 4]]) 

160 >>> B = matrix([[-2, 4], [5, 9]]) 

161 >>> A + B 

162 matrix( 

163 [['-1.0', '6.0'], 

164 ['8.0', '13.0']]) 

165 >>> A - B 

166 matrix( 

167 [['3.0', '-2.0'], 

168 ['-2.0', '-5.0']]) 

169 >>> A + ones(3) # doctest:+ELLIPSIS 

170 Traceback (most recent call last): 

171 ... 

172 ValueError: incompatible dimensions for addition 

173 

174 It is possible to multiply or add matrices and scalars. In the latter case the 

175 operation will be done element-wise:: 

176 

177 >>> A * 2 

178 matrix( 

179 [['2.0', '4.0'], 

180 ['6.0', '8.0']]) 

181 >>> A / 4 

182 matrix( 

183 [['0.25', '0.5'], 

184 ['0.75', '1.0']]) 

185 >>> A - 1 

186 matrix( 

187 [['0.0', '1.0'], 

188 ['2.0', '3.0']]) 

189 

190 Of course you can perform matrix multiplication, if the dimensions are 

191 compatible, using ``@`` (for Python >= 3.5) or ``*``. For clarity, ``@`` is 

192 recommended (`PEP 465 <https://www.python.org/dev/peps/pep-0465/>`), because 

193 the meaning of ``*`` is different in many other Python libraries such as NumPy. 

194 

195 >>> A @ B # doctest:+SKIP 

196 matrix( 

197 [['8.0', '22.0'], 

198 ['14.0', '48.0']]) 

199 >>> A * B # same as A @ B 

200 matrix( 

201 [['8.0', '22.0'], 

202 ['14.0', '48.0']]) 

203 >>> matrix([[1, 2, 3]]) * matrix([[-6], [7], [-2]]) 

204 matrix( 

205 [['2.0']]) 

206 

207 .. 

208 COMMENT: TODO: the above "doctest:+SKIP" may be removed as soon as we 

209 have dropped support for Python 3.4 and below. 

210 

211 You can raise powers of square matrices:: 

212 

213 >>> A**2 

214 matrix( 

215 [['7.0', '10.0'], 

216 ['15.0', '22.0']]) 

217 

218 Negative powers will calculate the inverse:: 

219 

220 >>> A**-1 

221 matrix( 

222 [['-2.0', '1.0'], 

223 ['1.5', '-0.5']]) 

224 >>> A * A**-1 

225 matrix( 

226 [['1.0', '1.0842021724855e-19'], 

227 ['-2.16840434497101e-19', '1.0']]) 

228 

229 

230 

231 Matrix transposition is straightforward:: 

232 

233 >>> A = ones(2, 3) 

234 >>> A 

235 matrix( 

236 [['1.0', '1.0', '1.0'], 

237 ['1.0', '1.0', '1.0']]) 

238 >>> A.T 

239 matrix( 

240 [['1.0', '1.0'], 

241 ['1.0', '1.0'], 

242 ['1.0', '1.0']]) 

243 

244 Norms 

245 ..... 

246 

247 Sometimes you need to know how "large" a matrix or vector is. Due to their 

248 multidimensional nature it's not possible to compare them, but there are 

249 several functions to map a matrix or a vector to a positive real number, the 

250 so called norms. 

251 

252 For vectors the p-norm is intended, usually the 1-, the 2- and the oo-norm are 

253 used. 

254 

255 >>> x = matrix([-10, 2, 100]) 

256 >>> norm(x, 1) 

257 mpf('112.0') 

258 >>> norm(x, 2) 

259 mpf('100.5186549850325') 

260 >>> norm(x, inf) 

261 mpf('100.0') 

262 

263 Please note that the 2-norm is the most used one, though it is more expensive 

264 to calculate than the 1- or oo-norm. 

265 

266 It is possible to generalize some vector norms to matrix norm:: 

267 

268 >>> A = matrix([[1, -1000], [100, 50]]) 

269 >>> mnorm(A, 1) 

270 mpf('1050.0') 

271 >>> mnorm(A, inf) 

272 mpf('1001.0') 

273 >>> mnorm(A, 'F') 

274 mpf('1006.2310867787777') 

275 

276 The last norm (the "Frobenius-norm") is an approximation for the 2-norm, which 

277 is hard to calculate and not available. The Frobenius-norm lacks some 

278 mathematical properties you might expect from a norm. 

279 """ 

280 

281 def __init__(self, *args, **kwargs): 

282 self.__data = {} 

283 # LU decompostion cache, this is useful when solving the same system 

284 # multiple times, when calculating the inverse and when calculating the 

285 # determinant 

286 self._LU = None 

287 if "force_type" in kwargs: 

288 warnings.warn("The force_type argument was removed, it did not work" 

289 " properly anyway. If you want to force floating-point or" 

290 " interval computations, use the respective methods from `fp`" 

291 " or `mp` instead, e.g., `fp.matrix()` or `iv.matrix()`." 

292 " If you want to truncate values to integer, use .apply(int) instead.") 

293 if isinstance(args[0], (list, tuple)): 

294 if isinstance(args[0][0], (list, tuple)): 

295 # interpret nested list as matrix 

296 A = args[0] 

297 self.__rows = len(A) 

298 self.__cols = len(A[0]) 

299 for i, row in enumerate(A): 

300 for j, a in enumerate(row): 

301 # note: this will call __setitem__ which will call self.ctx.convert() to convert the datatype. 

302 self[i, j] = a 

303 else: 

304 # interpret list as row vector 

305 v = args[0] 

306 self.__rows = len(v) 

307 self.__cols = 1 

308 for i, e in enumerate(v): 

309 self[i, 0] = e 

310 elif isinstance(args[0], int): 

311 # create empty matrix of given dimensions 

312 if len(args) == 1: 

313 self.__rows = self.__cols = args[0] 

314 else: 

315 if not isinstance(args[1], int): 

316 raise TypeError("expected int") 

317 self.__rows = args[0] 

318 self.__cols = args[1] 

319 elif isinstance(args[0], _matrix): 

320 A = args[0] 

321 self.__rows = A._matrix__rows 

322 self.__cols = A._matrix__cols 

323 for i in xrange(A.__rows): 

324 for j in xrange(A.__cols): 

325 self[i, j] = A[i, j] 

326 elif hasattr(args[0], 'tolist'): 

327 A = self.ctx.matrix(args[0].tolist()) 

328 self.__data = A._matrix__data 

329 self.__rows = A._matrix__rows 

330 self.__cols = A._matrix__cols 

331 else: 

332 raise TypeError('could not interpret given arguments') 

333 

334 def apply(self, f): 

335 """ 

336 Return a copy of self with the function `f` applied elementwise. 

337 """ 

338 new = self.ctx.matrix(self.__rows, self.__cols) 

339 for i in xrange(self.__rows): 

340 for j in xrange(self.__cols): 

341 new[i,j] = f(self[i,j]) 

342 return new 

343 

344 def __nstr__(self, n=None, **kwargs): 

345 # Build table of string representations of the elements 

346 res = [] 

347 # Track per-column max lengths for pretty alignment 

348 maxlen = [0] * self.cols 

349 for i in range(self.rows): 

350 res.append([]) 

351 for j in range(self.cols): 

352 if n: 

353 string = self.ctx.nstr(self[i,j], n, **kwargs) 

354 else: 

355 string = str(self[i,j]) 

356 res[-1].append(string) 

357 maxlen[j] = max(len(string), maxlen[j]) 

358 # Patch strings together 

359 for i, row in enumerate(res): 

360 for j, elem in enumerate(row): 

361 # Pad each element up to maxlen so the columns line up 

362 row[j] = elem.rjust(maxlen[j]) 

363 res[i] = "[" + colsep.join(row) + "]" 

364 return rowsep.join(res) 

365 

366 def __str__(self): 

367 return self.__nstr__() 

368 

369 def _toliststr(self, avoid_type=False): 

370 """ 

371 Create a list string from a matrix. 

372 

373 If avoid_type: avoid multiple 'mpf's. 

374 """ 

375 # XXX: should be something like self.ctx._types 

376 typ = self.ctx.mpf 

377 s = '[' 

378 for i in xrange(self.__rows): 

379 s += '[' 

380 for j in xrange(self.__cols): 

381 if not avoid_type or not isinstance(self[i,j], typ): 

382 a = repr(self[i,j]) 

383 else: 

384 a = "'" + str(self[i,j]) + "'" 

385 s += a + ', ' 

386 s = s[:-2] 

387 s += '],\n ' 

388 s = s[:-3] 

389 s += ']' 

390 return s 

391 

392 def tolist(self): 

393 """ 

394 Convert the matrix to a nested list. 

395 """ 

396 return [[self[i,j] for j in range(self.__cols)] for i in range(self.__rows)] 

397 

398 def __repr__(self): 

399 if self.ctx.pretty: 

400 return self.__str__() 

401 s = 'matrix(\n' 

402 s += self._toliststr(avoid_type=True) + ')' 

403 return s 

404 

405 def __get_element(self, key): 

406 ''' 

407 Fast extraction of the i,j element from the matrix 

408 This function is for private use only because is unsafe: 

409 1. Does not check on the value of key it expects key to be a integer tuple (i,j) 

410 2. Does not check bounds 

411 ''' 

412 if key in self.__data: 

413 return self.__data[key] 

414 else: 

415 return self.ctx.zero 

416 

417 def __set_element(self, key, value): 

418 ''' 

419 Fast assignment of the i,j element in the matrix 

420 This function is unsafe: 

421 1. Does not check on the value of key it expects key to be a integer tuple (i,j) 

422 2. Does not check bounds 

423 3. Does not check the value type 

424 4. Does not reset the LU cache 

425 ''' 

426 if value: # only store non-zeros 

427 self.__data[key] = value 

428 elif key in self.__data: 

429 del self.__data[key] 

430 

431 

432 def __getitem__(self, key): 

433 ''' 

434 Getitem function for mp matrix class with slice index enabled 

435 it allows the following assingments 

436 scalar to a slice of the matrix 

437 B = A[:,2:6] 

438 ''' 

439 # Convert vector to matrix indexing 

440 if isinstance(key, int) or isinstance(key,slice): 

441 # only sufficent for vectors 

442 if self.__rows == 1: 

443 key = (0, key) 

444 elif self.__cols == 1: 

445 key = (key, 0) 

446 else: 

447 raise IndexError('insufficient indices for matrix') 

448 

449 if isinstance(key[0],slice) or isinstance(key[1],slice): 

450 

451 #Rows 

452 if isinstance(key[0],slice): 

453 #Check bounds 

454 if (key[0].start is None or key[0].start >= 0) and \ 

455 (key[0].stop is None or key[0].stop <= self.__rows+1): 

456 # Generate indices 

457 rows = xrange(*key[0].indices(self.__rows)) 

458 else: 

459 raise IndexError('Row index out of bounds') 

460 else: 

461 # Single row 

462 rows = [key[0]] 

463 

464 # Columns 

465 if isinstance(key[1],slice): 

466 # Check bounds 

467 if (key[1].start is None or key[1].start >= 0) and \ 

468 (key[1].stop is None or key[1].stop <= self.__cols+1): 

469 # Generate indices 

470 columns = xrange(*key[1].indices(self.__cols)) 

471 else: 

472 raise IndexError('Column index out of bounds') 

473 

474 else: 

475 # Single column 

476 columns = [key[1]] 

477 

478 # Create matrix slice 

479 m = self.ctx.matrix(len(rows),len(columns)) 

480 

481 # Assign elements to the output matrix 

482 for i,x in enumerate(rows): 

483 for j,y in enumerate(columns): 

484 m.__set_element((i,j),self.__get_element((x,y))) 

485 

486 return m 

487 

488 else: 

489 # single element extraction 

490 if key[0] >= self.__rows or key[1] >= self.__cols: 

491 raise IndexError('matrix index out of range') 

492 if key in self.__data: 

493 return self.__data[key] 

494 else: 

495 return self.ctx.zero 

496 

497 def __setitem__(self, key, value): 

498 # setitem function for mp matrix class with slice index enabled 

499 # it allows the following assingments 

500 # scalar to a slice of the matrix 

501 # A[:,2:6] = 2.5 

502 # submatrix to matrix (the value matrix should be the same size as the slice size) 

503 # A[3,:] = B where A is n x m and B is n x 1 

504 # Convert vector to matrix indexing 

505 if isinstance(key, int) or isinstance(key,slice): 

506 # only sufficent for vectors 

507 if self.__rows == 1: 

508 key = (0, key) 

509 elif self.__cols == 1: 

510 key = (key, 0) 

511 else: 

512 raise IndexError('insufficient indices for matrix') 

513 # Slice indexing 

514 if isinstance(key[0],slice) or isinstance(key[1],slice): 

515 # Rows 

516 if isinstance(key[0],slice): 

517 # Check bounds 

518 if (key[0].start is None or key[0].start >= 0) and \ 

519 (key[0].stop is None or key[0].stop <= self.__rows+1): 

520 # generate row indices 

521 rows = xrange(*key[0].indices(self.__rows)) 

522 else: 

523 raise IndexError('Row index out of bounds') 

524 else: 

525 # Single row 

526 rows = [key[0]] 

527 # Columns 

528 if isinstance(key[1],slice): 

529 # Check bounds 

530 if (key[1].start is None or key[1].start >= 0) and \ 

531 (key[1].stop is None or key[1].stop <= self.__cols+1): 

532 # Generate column indices 

533 columns = xrange(*key[1].indices(self.__cols)) 

534 else: 

535 raise IndexError('Column index out of bounds') 

536 else: 

537 # Single column 

538 columns = [key[1]] 

539 # Assign slice with a scalar 

540 if isinstance(value,self.ctx.matrix): 

541 # Assign elements to matrix if input and output dimensions match 

542 if len(rows) == value.rows and len(columns) == value.cols: 

543 for i,x in enumerate(rows): 

544 for j,y in enumerate(columns): 

545 self.__set_element((x,y), value.__get_element((i,j))) 

546 else: 

547 raise ValueError('Dimensions do not match') 

548 else: 

549 # Assign slice with scalars 

550 value = self.ctx.convert(value) 

551 for i in rows: 

552 for j in columns: 

553 self.__set_element((i,j), value) 

554 else: 

555 # Single element assingment 

556 # Check bounds 

557 if key[0] >= self.__rows or key[1] >= self.__cols: 

558 raise IndexError('matrix index out of range') 

559 # Convert and store value 

560 value = self.ctx.convert(value) 

561 if value: # only store non-zeros 

562 self.__data[key] = value 

563 elif key in self.__data: 

564 del self.__data[key] 

565 

566 if self._LU: 

567 self._LU = None 

568 return 

569 

570 def __iter__(self): 

571 for i in xrange(self.__rows): 

572 for j in xrange(self.__cols): 

573 yield self[i,j] 

574 

575 def __mul__(self, other): 

576 if isinstance(other, self.ctx.matrix): 

577 # dot multiplication TODO: use Strassen's method? 

578 if self.__cols != other.__rows: 

579 raise ValueError('dimensions not compatible for multiplication') 

580 new = self.ctx.matrix(self.__rows, other.__cols) 

581 for i in xrange(self.__rows): 

582 for j in xrange(other.__cols): 

583 new[i, j] = self.ctx.fdot((self[i,k], other[k,j]) 

584 for k in xrange(other.__rows)) 

585 return new 

586 else: 

587 # try scalar multiplication 

588 new = self.ctx.matrix(self.__rows, self.__cols) 

589 for i in xrange(self.__rows): 

590 for j in xrange(self.__cols): 

591 new[i, j] = other * self[i, j] 

592 return new 

593 

594 def __matmul__(self, other): 

595 return self.__mul__(other) 

596 

597 def __rmul__(self, other): 

598 # assume other is scalar and thus commutative 

599 if isinstance(other, self.ctx.matrix): 

600 raise TypeError("other should not be type of ctx.matrix") 

601 return self.__mul__(other) 

602 

603 def __pow__(self, other): 

604 # avoid cyclic import problems 

605 #from linalg import inverse 

606 if not isinstance(other, int): 

607 raise ValueError('only integer exponents are supported') 

608 if not self.__rows == self.__cols: 

609 raise ValueError('only powers of square matrices are defined') 

610 n = other 

611 if n == 0: 

612 return self.ctx.eye(self.__rows) 

613 if n < 0: 

614 n = -n 

615 neg = True 

616 else: 

617 neg = False 

618 i = n 

619 y = 1 

620 z = self.copy() 

621 while i != 0: 

622 if i % 2 == 1: 

623 y = y * z 

624 z = z*z 

625 i = i // 2 

626 if neg: 

627 y = self.ctx.inverse(y) 

628 return y 

629 

630 def __div__(self, other): 

631 # assume other is scalar and do element-wise divison 

632 assert not isinstance(other, self.ctx.matrix) 

633 new = self.ctx.matrix(self.__rows, self.__cols) 

634 for i in xrange(self.__rows): 

635 for j in xrange(self.__cols): 

636 new[i,j] = self[i,j] / other 

637 return new 

638 

639 __truediv__ = __div__ 

640 

641 def __add__(self, other): 

642 if isinstance(other, self.ctx.matrix): 

643 if not (self.__rows == other.__rows and self.__cols == other.__cols): 

644 raise ValueError('incompatible dimensions for addition') 

645 new = self.ctx.matrix(self.__rows, self.__cols) 

646 for i in xrange(self.__rows): 

647 for j in xrange(self.__cols): 

648 new[i,j] = self[i,j] + other[i,j] 

649 return new 

650 else: 

651 # assume other is scalar and add element-wise 

652 new = self.ctx.matrix(self.__rows, self.__cols) 

653 for i in xrange(self.__rows): 

654 for j in xrange(self.__cols): 

655 new[i,j] += self[i,j] + other 

656 return new 

657 

658 def __radd__(self, other): 

659 return self.__add__(other) 

660 

661 def __sub__(self, other): 

662 if isinstance(other, self.ctx.matrix) and not (self.__rows == other.__rows 

663 and self.__cols == other.__cols): 

664 raise ValueError('incompatible dimensions for subtraction') 

665 return self.__add__(other * (-1)) 

666 

667 def __pos__(self): 

668 """ 

669 +M returns a copy of M, rounded to current working precision. 

670 """ 

671 return (+1) * self 

672 

673 def __neg__(self): 

674 return (-1) * self 

675 

676 def __rsub__(self, other): 

677 return -self + other 

678 

679 def __eq__(self, other): 

680 return self.__rows == other.__rows and self.__cols == other.__cols \ 

681 and self.__data == other.__data 

682 

683 def __len__(self): 

684 if self.rows == 1: 

685 return self.cols 

686 elif self.cols == 1: 

687 return self.rows 

688 else: 

689 return self.rows # do it like numpy 

690 

691 def __getrows(self): 

692 return self.__rows 

693 

694 def __setrows(self, value): 

695 for key in self.__data.copy(): 

696 if key[0] >= value: 

697 del self.__data[key] 

698 self.__rows = value 

699 

700 rows = property(__getrows, __setrows, doc='number of rows') 

701 

702 def __getcols(self): 

703 return self.__cols 

704 

705 def __setcols(self, value): 

706 for key in self.__data.copy(): 

707 if key[1] >= value: 

708 del self.__data[key] 

709 self.__cols = value 

710 

711 cols = property(__getcols, __setcols, doc='number of columns') 

712 

713 def transpose(self): 

714 new = self.ctx.matrix(self.__cols, self.__rows) 

715 for i in xrange(self.__rows): 

716 for j in xrange(self.__cols): 

717 new[j,i] = self[i,j] 

718 return new 

719 

720 T = property(transpose) 

721 

722 def conjugate(self): 

723 return self.apply(self.ctx.conj) 

724 

725 def transpose_conj(self): 

726 return self.conjugate().transpose() 

727 

728 H = property(transpose_conj) 

729 

730 def copy(self): 

731 new = self.ctx.matrix(self.__rows, self.__cols) 

732 new.__data = self.__data.copy() 

733 return new 

734 

735 __copy__ = copy 

736 

737 def column(self, n): 

738 m = self.ctx.matrix(self.rows, 1) 

739 for i in range(self.rows): 

740 m[i] = self[i,n] 

741 return m 

742 

743class MatrixMethods(object): 

744 

745 def __init__(ctx): 

746 # XXX: subclass 

747 ctx.matrix = type('matrix', (_matrix,), {}) 

748 ctx.matrix.ctx = ctx 

749 ctx.matrix.convert = ctx.convert 

750 

751 def eye(ctx, n, **kwargs): 

752 """ 

753 Create square identity matrix n x n. 

754 """ 

755 A = ctx.matrix(n, **kwargs) 

756 for i in xrange(n): 

757 A[i,i] = 1 

758 return A 

759 

760 def diag(ctx, diagonal, **kwargs): 

761 """ 

762 Create square diagonal matrix using given list. 

763 

764 Example: 

765 >>> from mpmath import diag, mp 

766 >>> mp.pretty = False 

767 >>> diag([1, 2, 3]) 

768 matrix( 

769 [['1.0', '0.0', '0.0'], 

770 ['0.0', '2.0', '0.0'], 

771 ['0.0', '0.0', '3.0']]) 

772 """ 

773 A = ctx.matrix(len(diagonal), **kwargs) 

774 for i in xrange(len(diagonal)): 

775 A[i,i] = diagonal[i] 

776 return A 

777 

778 def zeros(ctx, *args, **kwargs): 

779 """ 

780 Create matrix m x n filled with zeros. 

781 One given dimension will create square matrix n x n. 

782 

783 Example: 

784 >>> from mpmath import zeros, mp 

785 >>> mp.pretty = False 

786 >>> zeros(2) 

787 matrix( 

788 [['0.0', '0.0'], 

789 ['0.0', '0.0']]) 

790 """ 

791 if len(args) == 1: 

792 m = n = args[0] 

793 elif len(args) == 2: 

794 m = args[0] 

795 n = args[1] 

796 else: 

797 raise TypeError('zeros expected at most 2 arguments, got %i' % len(args)) 

798 A = ctx.matrix(m, n, **kwargs) 

799 for i in xrange(m): 

800 for j in xrange(n): 

801 A[i,j] = 0 

802 return A 

803 

804 def ones(ctx, *args, **kwargs): 

805 """ 

806 Create matrix m x n filled with ones. 

807 One given dimension will create square matrix n x n. 

808 

809 Example: 

810 >>> from mpmath import ones, mp 

811 >>> mp.pretty = False 

812 >>> ones(2) 

813 matrix( 

814 [['1.0', '1.0'], 

815 ['1.0', '1.0']]) 

816 """ 

817 if len(args) == 1: 

818 m = n = args[0] 

819 elif len(args) == 2: 

820 m = args[0] 

821 n = args[1] 

822 else: 

823 raise TypeError('ones expected at most 2 arguments, got %i' % len(args)) 

824 A = ctx.matrix(m, n, **kwargs) 

825 for i in xrange(m): 

826 for j in xrange(n): 

827 A[i,j] = 1 

828 return A 

829 

830 def hilbert(ctx, m, n=None): 

831 """ 

832 Create (pseudo) hilbert matrix m x n. 

833 One given dimension will create hilbert matrix n x n. 

834 

835 The matrix is very ill-conditioned and symmetric, positive definite if 

836 square. 

837 """ 

838 if n is None: 

839 n = m 

840 A = ctx.matrix(m, n) 

841 for i in xrange(m): 

842 for j in xrange(n): 

843 A[i,j] = ctx.one / (i + j + 1) 

844 return A 

845 

846 def randmatrix(ctx, m, n=None, min=0, max=1, **kwargs): 

847 """ 

848 Create a random m x n matrix. 

849 

850 All values are >= min and <max. 

851 n defaults to m. 

852 

853 Example: 

854 >>> from mpmath import randmatrix 

855 >>> randmatrix(2) # doctest:+SKIP 

856 matrix( 

857 [['0.53491598236191806', '0.57195669543302752'], 

858 ['0.85589992269513615', '0.82444367501382143']]) 

859 """ 

860 if not n: 

861 n = m 

862 A = ctx.matrix(m, n, **kwargs) 

863 for i in xrange(m): 

864 for j in xrange(n): 

865 A[i,j] = ctx.rand() * (max - min) + min 

866 return A 

867 

868 def swap_row(ctx, A, i, j): 

869 """ 

870 Swap row i with row j. 

871 """ 

872 if i == j: 

873 return