Coverage for /usr/lib/python3/dist-packages/sympy/polys/matrices/sdm.py: 19%

1"""

3Module for the SDM class.

5"""

7from operator import add, neg, pos, sub, mul

8from collections import defaultdict

10from sympy.utilities.iterables import _strongly_connected_components

12from .exceptions import DMBadInputError, DMDomainError, DMShapeError

14from .ddm import DDM

15from .lll import ddm_lll, ddm_lll_transform

16from sympy.polys.domains import QQ

19class SDM(dict):

20 r"""Sparse matrix based on polys domain elements

22 This is a dict subclass and is a wrapper for a dict of dicts that supports

23 basic matrix arithmetic +, -, *, **.

26 In order to create a new :py:class:`~.SDM`, a dict

27 of dicts mapping non-zero elements to their

28 corresponding row and column in the matrix is needed.

30 We also need to specify the shape and :py:class:`~.Domain`

31 of our :py:class:`~.SDM` object.

33 We declare a 2x2 :py:class:`~.SDM` matrix belonging

34 to QQ domain as shown below.

35 The 2x2 Matrix in the example is

37 .. math::

38 A = \left[\begin{array}{ccc}

39 0 & \frac{1}{2} \\

40 0 & 0 \end{array} \right]

43 >>> from sympy.polys.matrices.sdm import SDM

44 >>> from sympy import QQ

45 >>> elemsdict = {0:{1:QQ(1, 2)}}

46 >>> A = SDM(elemsdict, (2, 2), QQ)

47 >>> A

48 {0: {1: 1/2}}

50 We can manipulate :py:class:`~.SDM` the same way

51 as a Matrix class

53 >>> from sympy import ZZ

54 >>> A = SDM({0:{1: ZZ(2)}, 1:{0:ZZ(1)}}, (2, 2), ZZ)

55 >>> B = SDM({0:{0: ZZ(3)}, 1:{1:ZZ(4)}}, (2, 2), ZZ)

56 >>> A + B

57 {0: {0: 3, 1: 2}, 1: {0: 1, 1: 4}}

59 Multiplication

61 >>> A*B

62 {0: {1: 8}, 1: {0: 3}}

63 >>> A*ZZ(2)

64 {0: {1: 4}, 1: {0: 2}}

66 """

68 fmt = 'sparse'

70 def __init__(self, elemsdict, shape, domain):

71 super().__init__(elemsdict)

72 self.shape = self.rows, self.cols = m, n = shape

73 self.domain = domain

75 if not all(0 <= r < m for r in self):

76 raise DMBadInputError("Row out of range")

77 if not all(0 <= c < n for row in self.values() for c in row):

78 raise DMBadInputError("Column out of range")

80 def getitem(self, i, j):

81 try:

82 return self[i][j]

83 except KeyError:

84 m, n = self.shape

85 if -m <= i < m and -n <= j < n:

86 try:

87 return self[i % m][j % n]

88 except KeyError:

89 return self.domain.zero

90 else:

91 raise IndexError("index out of range")

93 def setitem(self, i, j, value):

94 m, n = self.shape

95 if not (-m <= i < m and -n <= j < n):

96 raise IndexError("index out of range")

97 i, j = i % m, j % n

98 if value:

99 try:

100 self[i][j] = value

101 except KeyError:

102 self[i] = {j: value}

103 else:

104 rowi = self.get(i, None)

105 if rowi is not None:

106 try:

107 del rowi[j]

108 except KeyError:

109 pass

110 else:

111 if not rowi:

112 del self[i]

113

114 def extract_slice(self, slice1, slice2):

115 m, n = self.shape

116 ri = range(m)[slice1]

117 ci = range(n)[slice2]

118

119 sdm = {}

120 for i, row in self.items():

121 if i in ri:

122 row = {ci.index(j): e for j, e in row.items() if j in ci}

123 if row:

124 sdm[ri.index(i)] = row

125

126 return self.new(sdm, (len(ri), len(ci)), self.domain)

127

128 def extract(self, rows, cols):

129 if not (self and rows and cols):

130 return self.zeros((len(rows), len(cols)), self.domain)

131

132 m, n = self.shape

133 if not (-m <= min(rows) <= max(rows) < m):

134 raise IndexError('Row index out of range')

135 if not (-n <= min(cols) <= max(cols) < n):

136 raise IndexError('Column index out of range')

137

138 # rows and cols can contain duplicates e.g. M[[1, 2, 2], [0, 1]]

139 # Build a map from row/col in self to list of rows/cols in output

140 rowmap = defaultdict(list)

141 colmap = defaultdict(list)

142 for i2, i1 in enumerate(rows):

143 rowmap[i1 % m].append(i2)

144 for j2, j1 in enumerate(cols):

145 colmap[j1 % n].append(j2)

146

147 # Used to efficiently skip zero rows/cols

148 rowset = set(rowmap)

149 colset = set(colmap)

150

151 sdm1 = self

152 sdm2 = {}

153 for i1 in rowset & set(sdm1):

154 row1 = sdm1[i1]

155 row2 = {}

156 for j1 in colset & set(row1):

157 row1_j1 = row1[j1]

158 for j2 in colmap[j1]:

159 row2[j2] = row1_j1

160 if row2:

161 for i2 in rowmap[i1]:

162 sdm2[i2] = row2.copy()

163

164 return self.new(sdm2, (len(rows), len(cols)), self.domain)

165

166 def __str__(self):

167 rowsstr = []

168 for i, row in self.items():

169 elemsstr = ', '.join('%s: %s' % (j, elem) for j, elem in row.items())

170 rowsstr.append('%s: {%s}' % (i, elemsstr))

171 return '{%s}' % ', '.join(rowsstr)

172

173 def __repr__(self):

174 cls = type(self).__name__

175 rows = dict.__repr__(self)

176 return '%s(%s, %s, %s)' % (cls, rows, self.shape, self.domain)

177

178 @classmethod

179 def new(cls, sdm, shape, domain):

180 """

181

182 Parameters

183 ==========

184

185 sdm: A dict of dicts for non-zero elements in SDM

186 shape: tuple representing dimension of SDM

187 domain: Represents :py:class:`~.Domain` of SDM

188

189 Returns

190 =======

191

192 An :py:class:`~.SDM` object

193

194 Examples

195 ========

196

197 >>> from sympy.polys.matrices.sdm import SDM

198 >>> from sympy import QQ

199 >>> elemsdict = {0:{1: QQ(2)}}

200 >>> A = SDM.new(elemsdict, (2, 2), QQ)

201 >>> A

202 {0: {1: 2}}

203

204 """

205 return cls(sdm, shape, domain)

206

207 def copy(A):

208 """

209 Returns the copy of a :py:class:`~.SDM` object

210

211 Examples

212 ========

213

214 >>> from sympy.polys.matrices.sdm import SDM

215 >>> from sympy import QQ

216 >>> elemsdict = {0:{1:QQ(2)}, 1:{}}

217 >>> A = SDM(elemsdict, (2, 2), QQ)

218 >>> B = A.copy()

219 >>> B

220 {0: {1: 2}, 1: {}}

221

222 """

223 Ac = {i: Ai.copy() for i, Ai in A.items()}

224 return A.new(Ac, A.shape, A.domain)

225

226 @classmethod

227 def from_list(cls, ddm, shape, domain):

228 """

229

230 Parameters

231 ==========

232

233 ddm:

234 list of lists containing domain elements

235 shape:

236 Dimensions of :py:class:`~.SDM` matrix

237 domain:

238 Represents :py:class:`~.Domain` of :py:class:`~.SDM` object

239

240 Returns

241 =======

242

243 :py:class:`~.SDM` containing elements of ddm

244

245 Examples

246 ========

247

248 >>> from sympy.polys.matrices.sdm import SDM

249 >>> from sympy import QQ

250 >>> ddm = [[QQ(1, 2), QQ(0)], [QQ(0), QQ(3, 4)]]

251 >>> A = SDM.from_list(ddm, (2, 2), QQ)

252 >>> A

253 {0: {0: 1/2}, 1: {1: 3/4}}

254

255 """

256

257 m, n = shape

258 if not (len(ddm) == m and all(len(row) == n for row in ddm)):

259 raise DMBadInputError("Inconsistent row-list/shape")

260 getrow = lambda i: {j:ddm[i][j] for j in range(n) if ddm[i][j]}

261 irows = ((i, getrow(i)) for i in range(m))

262 sdm = {i: row for i, row in irows if row}

263 return cls(sdm, shape, domain)

264

265 @classmethod

266 def from_ddm(cls, ddm):

267 """

268 converts object of :py:class:`~.DDM` to

269 :py:class:`~.SDM`

270

271 Examples

272 ========

273

274 >>> from sympy.polys.matrices.ddm import DDM

275 >>> from sympy.polys.matrices.sdm import SDM

276 >>> from sympy import QQ

277 >>> ddm = DDM( [[QQ(1, 2), 0], [0, QQ(3, 4)]], (2, 2), QQ)

278 >>> A = SDM.from_ddm(ddm)

279 >>> A

280 {0: {0: 1/2}, 1: {1: 3/4}}

281

282 """

283 return cls.from_list(ddm, ddm.shape, ddm.domain)

284

285 def to_list(M):

286 """

287

288 Converts a :py:class:`~.SDM` object to a list

289

290 Examples

291 ========

292

293 >>> from sympy.polys.matrices.sdm import SDM

294 >>> from sympy import QQ

295 >>> elemsdict = {0:{1:QQ(2)}, 1:{}}

296 >>> A = SDM(elemsdict, (2, 2), QQ)

297 >>> A.to_list()

298 [[0, 2], [0, 0]]

299

300 """

301 m, n = M.shape

302 zero = M.domain.zero

303 ddm = [[zero] * n for _ in range(m)]

304 for i, row in M.items():

305 for j, e in row.items():

306 ddm[i][j] = e

307 return ddm

308

309 def to_list_flat(M):

310 m, n = M.shape

311 zero = M.domain.zero

312 flat = [zero] * (m * n)

313 for i, row in M.items():

314 for j, e in row.items():

315 flat[i*n + j] = e

316 return flat

317

318 def to_dok(M):

319 return {(i, j): e for i, row in M.items() for j, e in row.items()}

320

321 def to_ddm(M):

322 """

323 Convert a :py:class:`~.SDM` object to a :py:class:`~.DDM` object

324

325 Examples

326 ========

327

328 >>> from sympy.polys.matrices.sdm import SDM

329 >>> from sympy import QQ

330 >>> A = SDM({0:{1:QQ(2)}, 1:{}}, (2, 2), QQ)

331 >>> A.to_ddm()

332 [[0, 2], [0, 0]]

333

334 """

335 return DDM(M.to_list(), M.shape, M.domain)

336

337 def to_sdm(M):

338 return M

339

340 @classmethod

341 def zeros(cls, shape, domain):

342 r"""

343

344 Returns a :py:class:`~.SDM` of size shape,

345 belonging to the specified domain

346

347 In the example below we declare a matrix A where,

348

349 .. math::

350 A := \left[\begin{array}{ccc}

351 0 & 0 & 0 \\

352 0 & 0 & 0 \end{array} \right]

353

354 >>> from sympy.polys.matrices.sdm import SDM

355 >>> from sympy import QQ

356 >>> A = SDM.zeros((2, 3), QQ)

357 >>> A

358 {}

359

360 """

361 return cls({}, shape, domain)

362

363 @classmethod

364 def ones(cls, shape, domain):

365 one = domain.one

366 m, n = shape

367 row = dict(zip(range(n), [one]*n))

368 sdm = {i: row.copy() for i in range(m)}

369 return cls(sdm, shape, domain)

370

371 @classmethod

372 def eye(cls, shape, domain):

373 """

374

375 Returns a identity :py:class:`~.SDM` matrix of dimensions

376 size x size, belonging to the specified domain

377

378 Examples

379 ========

380

381 >>> from sympy.polys.matrices.sdm import SDM

382 >>> from sympy import QQ

383 >>> I = SDM.eye((2, 2), QQ)

384 >>> I

385 {0: {0: 1}, 1: {1: 1}}

386

387 """

388 rows, cols = shape

389 one = domain.one

390 sdm = {i: {i: one} for i in range(min(rows, cols))}

391 return cls(sdm, shape, domain)

392

393 @classmethod

394 def diag(cls, diagonal, domain, shape):

395 sdm = {i: {i: v} for i, v in enumerate(diagonal) if v}

396 return cls(sdm, shape, domain)

397

398 def transpose(M):

399 """

400

401 Returns the transpose of a :py:class:`~.SDM` matrix

402

403 Examples

404 ========

405

406 >>> from sympy.polys.matrices.sdm import SDM

407 >>> from sympy import QQ

408 >>> A = SDM({0:{1:QQ(2)}, 1:{}}, (2, 2), QQ)

409 >>> A.transpose()

410 {1: {0: 2}}

411

412 """

413 MT = sdm_transpose(M)

414 return M.new(MT, M.shape[::-1], M.domain)

415

416 def __add__(A, B):

417 if not isinstance(B, SDM):

418 return NotImplemented

419 return A.add(B)

420

421 def __sub__(A, B):

422 if not isinstance(B, SDM):

423 return NotImplemented

424 return A.sub(B)

425

426 def __neg__(A):

427 return A.neg()

428

429 def __mul__(A, B):

430 """A * B"""

431 if isinstance(B, SDM):

432 return A.matmul(B)

433 elif B in A.domain:

434 return A.mul(B)

435 else:

436 return NotImplemented

437

438 def __rmul__(a, b):

439 if b in a.domain:

440 return a.rmul(b)

441 else:

442 return NotImplemented

443

444 def matmul(A, B):

445 """

446 Performs matrix multiplication of two SDM matrices

447

448 Parameters

449 ==========

450

451 A, B: SDM to multiply

452

453 Returns

454 =======

455

456 SDM

457 SDM after multiplication

458

459 Raises

460 ======

461

462 DomainError

463 If domain of A does not match

464 with that of B

465

466 Examples

467 ========

468

469 >>> from sympy import ZZ

470 >>> from sympy.polys.matrices.sdm import SDM

471 >>> A = SDM({0:{1: ZZ(2)}, 1:{0:ZZ(1)}}, (2, 2), ZZ)

472 >>> B = SDM({0:{0:ZZ(2), 1:ZZ(3)}, 1:{0:ZZ(4)}}, (2, 2), ZZ)

473 >>> A.matmul(B)

474 {0: {0: 8}, 1: {0: 2, 1: 3}}

475

476 """

477 if A.domain != B.domain:

478 raise DMDomainError

479 m, n = A.shape

480 n2, o = B.shape

481 if n != n2:

482 raise DMShapeError

483 C = sdm_matmul(A, B, A.domain, m, o)

484 return A.new(C, (m, o), A.domain)

485

486 def mul(A, b):

487 """

488 Multiplies each element of A with a scalar b

489

490 Examples

491 ========

492

493 >>> from sympy import ZZ

494 >>> from sympy.polys.matrices.sdm import SDM

495 >>> A = SDM({0:{1: ZZ(2)}, 1:{0:ZZ(1)}}, (2, 2), ZZ)

496 >>> A.mul(ZZ(3))

497 {0: {1: 6}, 1: {0: 3}}

498

499 """

500 Csdm = unop_dict(A, lambda aij: aij*b)

501 return A.new(Csdm, A.shape, A.domain)

502

503 def rmul(A, b):

504 Csdm = unop_dict(A, lambda aij: b*aij)

505 return A.new(Csdm, A.shape, A.domain)

506

507 def mul_elementwise(A, B):

508 if A.domain != B.domain:

509 raise DMDomainError

510 if A.shape != B.shape:

511 raise DMShapeError

512 zero = A.domain.zero

513 fzero = lambda e: zero

514 Csdm = binop_dict(A, B, mul, fzero, fzero)

515 return A.new(Csdm, A.shape, A.domain)

516

517 def add(A, B):

518 """

519

520 Adds two :py:class:`~.SDM` matrices

521

522 Examples

523 ========

524

525 >>> from sympy import ZZ

526 >>> from sympy.polys.matrices.sdm import SDM

527 >>> A = SDM({0:{1: ZZ(2)}, 1:{0:ZZ(1)}}, (2, 2), ZZ)

528 >>> B = SDM({0:{0: ZZ(3)}, 1:{1:ZZ(4)}}, (2, 2), ZZ)

529 >>> A.add(B)

530 {0: {0: 3, 1: 2}, 1: {0: 1, 1: 4}}

531

532 """

533

534 Csdm = binop_dict(A, B, add, pos, pos)

535 return A.new(Csdm, A.shape, A.domain)

536

537 def sub(A, B):

538 """

539

540 Subtracts two :py:class:`~.SDM` matrices

541

542 Examples

543 ========

544

545 >>> from sympy import ZZ

546 >>> from sympy.polys.matrices.sdm import SDM

547 >>> A = SDM({0:{1: ZZ(2)}, 1:{0:ZZ(1)}}, (2, 2), ZZ)

548 >>> B = SDM({0:{0: ZZ(3)}, 1:{1:ZZ(4)}}, (2, 2), ZZ)

549 >>> A.sub(B)

550 {0: {0: -3, 1: 2}, 1: {0: 1, 1: -4}}

551

552 """

553 Csdm = binop_dict(A, B, sub, pos, neg)

554 return A.new(Csdm, A.shape, A.domain)

555

556 def neg(A):

557 """

558

559 Returns the negative of a :py:class:`~.SDM` matrix

560

561 Examples

562 ========

563

564 >>> from sympy import ZZ

565 >>> from sympy.polys.matrices.sdm import SDM

566 >>> A = SDM({0:{1: ZZ(2)}, 1:{0:ZZ(1)}}, (2, 2), ZZ)

567 >>> A.neg()

568 {0: {1: -2}, 1: {0: -1}}

569

570 """

571 Csdm = unop_dict(A, neg)

572 return A.new(Csdm, A.shape, A.domain)

573

574 def convert_to(A, K):

575 """

576

577 Converts the :py:class:`~.Domain` of a :py:class:`~.SDM` matrix to K

578

579 Examples

580 ========

581

582 >>> from sympy import ZZ, QQ

583 >>> from sympy.polys.matrices.sdm import SDM

584 >>> A = SDM({0:{1: ZZ(2)}, 1:{0:ZZ(1)}}, (2, 2), ZZ)

585 >>> A.convert_to(QQ)

586 {0: {1: 2}, 1: {0: 1}}

587

588 """

589 Kold = A.domain

590 if K == Kold:

591 return A.copy()

592 Ak = unop_dict(A, lambda e: K.convert_from(e, Kold))

593 return A.new(Ak, A.shape, K)

594

595 def scc(A):

596 """Strongly connected components of a square matrix *A*.

597

598 Examples

599 ========

600

601 >>> from sympy import ZZ

602 >>> from sympy.polys.matrices.sdm import SDM

603 >>> A = SDM({0:{0: ZZ(2)}, 1:{1:ZZ(1)}}, (2, 2), ZZ)

604 >>> A.scc()

605 [[0], [1]]

606

607 See also

608 ========

609

610 sympy.polys.matrices.domainmatrix.DomainMatrix.scc

611 """

612 rows, cols = A.shape

613 assert rows == cols

614 V = range(rows)

615 Emap = {v: list(A.get(v, [])) for v in V}

616 return _strongly_connected_components(V, Emap)

617

618 def rref(A):

619 """

620

621 Returns reduced-row echelon form and list of pivots for the :py:class:`~.SDM`

622

623 Examples

624 ========

625

626 >>> from sympy import QQ

627 >>> from sympy.polys.matrices.sdm import SDM

628 >>> A = SDM({0:{0:QQ(1), 1:QQ(2)}, 1:{0:QQ(2), 1:QQ(4)}}, (2, 2), QQ)

629 >>> A.rref()

630 ({0: {0: 1, 1: 2}}, [0])

631

632 """

633 B, pivots, _ = sdm_irref(A)

634 return A.new(B, A.shape, A.domain), pivots

635

636 def inv(A):

637 """

638

639 Returns inverse of a matrix A

640

641 Examples

642 ========

643

644 >>> from sympy import QQ

645 >>> from sympy.polys.matrices.sdm import SDM

646 >>> A = SDM({0:{0:QQ(1), 1:QQ(2)}, 1:{0:QQ(3), 1:QQ(4)}}, (2, 2), QQ)

647 >>> A.inv()

648 {0: {0: -2, 1: 1}, 1: {0: 3/2, 1: -1/2}}

649

650 """

651 return A.from_ddm(A.to_ddm().inv())

652

653 def det(A):

654 """

655 Returns determinant of A

656

657 Examples

658 ========

659

660 >>> from sympy import QQ

661 >>> from sympy.polys.matrices.sdm import SDM

662 >>> A = SDM({0:{0:QQ(1), 1:QQ(2)}, 1:{0:QQ(3), 1:QQ(4)}}, (2, 2), QQ)

663 >>> A.det()

664 -2

665

666 """

667 return A.to_ddm().det()

668

669 def lu(A):

670 """

671

672 Returns LU decomposition for a matrix A

673

674 Examples

675 ========

676

677 >>> from sympy import QQ

678 >>> from sympy.polys.matrices.sdm import SDM

679 >>> A = SDM({0:{0:QQ(1), 1:QQ(2)}, 1:{0:QQ(3), 1:QQ(4)}}, (2, 2), QQ)

680 >>> A.lu()

681 ({0: {0: 1}, 1: {0: 3, 1: 1}}, {0: {0: 1, 1: 2}, 1: {1: -2}}, [])

682

683 """

684 L, U, swaps = A.to_ddm().lu()

685 return A.from_ddm(L), A.from_ddm(U), swaps

686

687 def lu_solve(A, b):

688 """

689

690 Uses LU decomposition to solve Ax = b,

691

692 Examples

693 ========

694

695 >>> from sympy import QQ

696 >>> from sympy.polys.matrices.sdm import SDM

697 >>> A = SDM({0:{0:QQ(1), 1:QQ(2)}, 1:{0:QQ(3), 1:QQ(4)}}, (2, 2), QQ)

698 >>> b = SDM({0:{0:QQ(1)}, 1:{0:QQ(2)}}, (2, 1), QQ)

699 >>> A.lu_solve(b)

700 {1: {0: 1/2}}

701

702 """

703 return A.from_ddm(A.to_ddm().lu_solve(b.to_ddm()))

704

705 def nullspace(A):

706 """

707

708 Returns nullspace for a :py:class:`~.SDM` matrix A

709

710 Examples

711 ========

712

713 >>> from sympy import QQ

714 >>> from sympy.polys.matrices.sdm import SDM

715 >>> A = SDM({0:{0:QQ(1), 1:QQ(2)}, 1:{0: QQ(2), 1: QQ(4)}}, (2, 2), QQ)

716 >>> A.nullspace()

717 ({0: {0: -2, 1: 1}}, [1])

718

719 """

720 ncols = A.shape[1]

721 one = A.domain.one

722 B, pivots, nzcols = sdm_irref(A)

723 K, nonpivots = sdm_nullspace_from_rref(B, one, ncols, pivots, nzcols)

724 K = dict(enumerate(K))

725 shape = (len(K), ncols)

726 return A.new(K, shape, A.domain), nonpivots

727

728 def particular(A):

729 ncols = A.shape[1]

730 B, pivots, nzcols = sdm_irref(A)

731 P = sdm_particular_from_rref(B, ncols, pivots)

732 rep = {0:P} if P else {}

733 return A.new(rep, (1, ncols-1), A.domain)

734

735 def hstack(A, *B):

736 """Horizontally stacks :py:class:`~.SDM` matrices.

737

738 Examples

739 ========

740

741 >>> from sympy import ZZ

742 >>> from sympy.polys.matrices.sdm import SDM

743

744 >>> A = SDM({0: {0: ZZ(1), 1: ZZ(2)}, 1: {0: ZZ(3), 1: ZZ(4)}}, (2, 2), ZZ)

745 >>> B = SDM({0: {0: ZZ(5), 1: ZZ(6)}, 1: {0: ZZ(7), 1: ZZ(8)}}, (2, 2), ZZ)

746 >>> A.hstack(B)

747 {0: {0: 1, 1: 2, 2: 5, 3: 6}, 1: {0: 3, 1: 4, 2: 7, 3: 8}}

748

749 >>> C = SDM({0: {0: ZZ(9), 1: ZZ(10)}, 1: {0: ZZ(11), 1: ZZ(12)}}, (2, 2), ZZ)

750 >>> A.hstack(B, C)

751 {0: {0: 1, 1: 2, 2: 5, 3: 6, 4: 9, 5: 10}, 1: {0: 3, 1: 4, 2: 7, 3: 8, 4: 11, 5: 12}}

752 """

753 Anew = dict(A.copy())

754 rows, cols = A.shape

755 domain = A.domain

756

757 for Bk in B:

758 Bkrows, Bkcols = Bk.shape

759 assert Bkrows == rows

760 assert Bk.domain == domain

761

762 for i, Bki in Bk.items():

763 Ai = Anew.get(i, None)

764 if Ai is None:

765 Anew[i] = Ai = {}

766 for j, Bkij in Bki.items():

767 Ai[j + cols] = Bkij

768 cols += Bkcols

769

770 return A.new(Anew, (rows, cols), A.domain)

771

772 def vstack(A, *B):

773 """Vertically stacks :py:class:`~.SDM` matrices.

774

775 Examples

776 ========

777

778 >>> from sympy import ZZ

779 >>> from sympy.polys.matrices.sdm import SDM

780

781 >>> A = SDM({0: {0: ZZ(1), 1: ZZ(2)}, 1: {0: ZZ(3), 1: ZZ(4)}}, (2, 2), ZZ)

782 >>> B = SDM({0: {0: ZZ(5), 1: ZZ(6)}, 1: {0: ZZ(7), 1: ZZ(8)}}, (2, 2), ZZ)

783 >>> A.vstack(B)

784 {0: {0: 1, 1: 2}, 1: {0: 3, 1: 4}, 2: {0: 5, 1: 6}, 3: {0: 7, 1: 8}}

785

786 >>> C = SDM({0: {0: ZZ(9), 1: ZZ(10)}, 1: {0: ZZ(11), 1: ZZ(12)}}, (2, 2), ZZ)

787 >>> A.vstack(B, C)

788 {0: {0: 1, 1: 2}, 1: {0: 3, 1: 4}, 2: {0: 5, 1: 6}, 3: {0: 7, 1: 8}, 4: {0: 9, 1: 10}, 5: {0: 11, 1: 12}}

789 """

790 Anew = dict(A.copy())

791 rows, cols = A.shape

792 domain = A.domain

793

794 for Bk in B:

795 Bkrows, Bkcols = Bk.shape

796 assert Bkcols == cols

797 assert Bk.domain == domain

798

799 for i, Bki in Bk.items():

800 Anew[i + rows] = Bki

801 rows += Bkrows

802

803 return A.new(Anew, (rows, cols), A.domain)

804

805 def applyfunc(self, func, domain):

806 sdm = {i: {j: func(e) for j, e in row.items()} for i, row in self.items()}

807 return self.new(sdm, self.shape, domain)

808

809 def charpoly(A):

810 """

811 Returns the coefficients of the characteristic polynomial

812 of the :py:class:`~.SDM` matrix. These elements will be domain elements.

813 The domain of the elements will be same as domain of the :py:class:`~.SDM`.

814

815 Examples

816 ========

817

818 >>> from sympy import QQ, Symbol

819 >>> from sympy.polys.matrices.sdm import SDM

820 >>> from sympy.polys import Poly

821 >>> A = SDM({0:{0:QQ(1), 1:QQ(2)}, 1:{0:QQ(3), 1:QQ(4)}}, (2, 2), QQ)

822 >>> A.charpoly()

823 [1, -5, -2]

824

825 We can create a polynomial using the

826 coefficients using :py:class:`~.Poly`

827

828 >>> x = Symbol('x')

829 >>> p = Poly(A.charpoly(), x, domain=A.domain)

830 >>> p

831 Poly(x**2 - 5*x - 2, x, domain='QQ')

832

833 """

834 return A.to_ddm().charpoly()

835

836 def is_zero_matrix(self):

837 """

838 Says whether this matrix has all zero entries.

839 """

840 return not self

841

842 def is_upper(self):

843 """

844 Says whether this matrix is upper-triangular. True can be returned

845 even if the matrix is not square.

846 """

847 return all(i <= j for i, row in self.items() for j in row)

848

849 def is_lower(self):

850 """

851 Says whether this matrix is lower-triangular. True can be returned

852 even if the matrix is not square.

853 """

854 return all(i >= j for i, row in self.items() for j in row)

855

856 def lll(A, delta=QQ(3, 4)):

857 return A.from_ddm(ddm_lll(A.to_ddm(), delta=delta))

858

859 def lll_transform(A, delta=QQ(3, 4)):

860 reduced, transform = ddm_lll_transform(A.to_ddm(), delta=delta)

861 return A.from_ddm(reduced), A.from_ddm(transform)

862

863

864def binop_dict(A, B, fab, fa, fb):

865 Anz, Bnz = set(A), set(B)

866 C = {}

867

868 for i in Anz & Bnz:

869 Ai, Bi = A[i], B[i]

870 Ci = {}

871 Anzi, Bnzi = set(Ai), set(Bi)

872 for j in Anzi & Bnzi:

873 Cij = fab(Ai[j], Bi[j])

874 if Cij:

875 Ci[j] = Cij

876 for j in Anzi - Bnzi:

877 Cij = fa(Ai[j])

878 if Cij:

879 Ci[j] = Cij

880 for j in Bnzi - Anzi:

881 Cij = fb(Bi[j])

882 if Cij:

883 Ci[j] = Cij

884 if Ci:

885 C[i] = Ci

886

887 for i in Anz - Bnz:

888 Ai = A[i]

889 Ci = {}

890 for j, Aij in Ai.items():

891 Cij = fa(Aij)

892 if Cij:

893 Ci[j] = Cij

894 if Ci:

895 C[i] = Ci

896

897 for i in Bnz - Anz:

898 Bi = B[i]

899 Ci = {}

900 for j, Bij in Bi.items():

901 Cij = fb(Bij)

902 if Cij:

903 Ci[j] = Cij

904 if Ci:

905 C[i] = Ci

906

907 return C

908

909

910def unop_dict(A, f):

911 B = {}

912 for i, Ai in A.items():

913 Bi = {}

914 for j, Aij in Ai.items():

915 Bij = f(Aij)

916 if Bij:

917 Bi[j] = Bij

918 if Bi:

919 B[i] = Bi

920 return B

921

922

923def sdm_transpose(M):

924 MT = {}

925 for i, Mi in M.items():

926 for j, Mij in Mi.items():

927 try:

928 MT[j][i] = Mij

929 except KeyError:

930 MT[j] = {i: Mij}

931 return MT

932

933

934def sdm_matmul(A, B, K, m, o):

935 #

936 # Should be fast if A and B are very sparse.

937 # Consider e.g. A = B = eye(1000).

938 #

939 # The idea here is that we compute C = A*B in terms of the rows of C and

940 # B since the dict of dicts representation naturally stores the matrix as

941 # rows. The ith row of C (Ci) is equal to the sum of Aik * Bk where Bk is

942 # the kth row of B. The algorithm below loops over each nonzero element

943 # Aik of A and if the corresponding row Bj is nonzero then we do

944 # Ci += Aik * Bk.

945 # To make this more efficient we don't need to loop over all elements Aik.

946 # Instead for each row Ai we compute the intersection of the nonzero

947 # columns in Ai with the nonzero rows in B. That gives the k such that

948 # Aik and Bk are both nonzero. In Python the intersection of two sets

949 # of int can be computed very efficiently.

950 #

951 if K.is_EXRAW:

952 return sdm_matmul_exraw(A, B, K, m, o)

953

954 C = {}

955 B_knz = set(B)

956 for i, Ai in A.items():

957 Ci = {}

958 Ai_knz = set(Ai)

959 for k in Ai_knz & B_knz:

960 Aik = Ai[k]

961 for j, Bkj in B[k].items():

962 Cij = Ci.get(j, None)

963 if Cij is not None:

964 Cij = Cij + Aik * Bkj

965 if Cij:

966 Ci[j] = Cij

967 else:

968 Ci.pop(j)

969 else:

970 Cij = Aik * Bkj

971 if Cij:

972 Ci[j] = Cij

973 if Ci:

974 C[i] = Ci

975 return C

976

977

978def sdm_matmul_exraw(A, B, K, m, o):

979 #

980 # Like sdm_matmul above except that:

981 #

982 # - Handles cases like 0*oo -> nan (sdm_matmul skips multipication by zero)

983 # - Uses K.sum (Add(*items)) for efficient addition of Expr

984 #

985 zero = K.zero

986 C = {}

987 B_knz = set(B)

988 for i, Ai in A.items():

989 Ci_list = defaultdict(list)

990 Ai_knz = set(Ai)

991

992 # Nonzero row/column pair

993 for k in Ai_knz & B_knz:

994 Aik = Ai[k]

995 if zero * Aik == zero:

996 # This is the main inner loop:

997 for j, Bkj in B[k].items():

998 Ci_list[j].append(Aik * Bkj)

999 else:

1000 for j in range(o):

1001 Ci_list[j].append(Aik * B[k].get(j, zero))

1002

1003 # Zero row in B, check for infinities in A

1004 for k in Ai_knz - B_knz:

1005 zAik = zero * Ai[k]

1006 if zAik != zero:

1007 for j in range(o):

1008 Ci_list[j].append(zAik)

1009

1010 # Add terms using K.sum (Add(*terms)) for efficiency

1011 Ci = {}

1012 for j, Cij_list in Ci_list.items():

1013 Cij = K.sum(Cij_list)

1014 if Cij:

1015 Ci[j] = Cij

1016 if Ci:

1017 C[i] = Ci

1018

1019 # Find all infinities in B

1020 for k, Bk in B.items():

1021 for j, Bkj in Bk.items():

1022 if zero * Bkj != zero:

1023 for i in range(m):

1024 Aik = A.get(i, {}).get(k, zero)

1025 # If Aik is not zero then this was handled above

1026 if Aik == zero:

1027 Ci = C.get(i, {})

1028 Cij = Ci.get(j, zero) + Aik * Bkj

1029 if Cij != zero:

1030 Ci[j] = Cij

1031 else: # pragma: no cover

1032 # Not sure how we could get here but let's raise an

1033 # exception just in case.

1034 raise RuntimeError

1035 C[i] = Ci

1036

1037 return C

1038

1039

1040def sdm_irref(A):

1041 """RREF and pivots of a sparse matrix *A*.

1042

1043 Compute the reduced row echelon form (RREF) of the matrix *A* and return a

1044 list of the pivot columns. This routine does not work in place and leaves

1045 the original matrix *A* unmodified.

1046

1047 Examples

1048 ========

1049

1050 This routine works with a dict of dicts sparse representation of a matrix:

1051

1052 >>> from sympy import QQ

1053 >>> from sympy.polys.matrices.sdm import sdm_irref

1054 >>> A = {0: {0: QQ(1), 1: QQ(2)}, 1: {0: QQ(3), 1: QQ(4)}}

1055 >>> Arref, pivots, _ = sdm_irref(A)

1056 >>> Arref

1057 {0: {0: 1}, 1: {1: 1}}

1058 >>> pivots

1059 [0, 1]

1060

1061 The analogous calculation with :py:class:`~.Matrix` would be

1062

1063 >>> from sympy import Matrix

1064 >>> M = Matrix([[1, 2], [3, 4]])

1065 >>> Mrref, pivots = M.rref()

1066 >>> Mrref

1067 Matrix([

1068 [1, 0],

1069 [0, 1]])

1070 >>> pivots

1071 (0, 1)

1072

1073 Notes

1074 =====

1075

1076 The cost of this algorithm is determined purely by the nonzero elements of

1077 the matrix. No part of the cost of any step in this algorithm depends on

1078 the number of rows or columns in the matrix. No step depends even on the

1079 number of nonzero rows apart from the primary loop over those rows. The

1080 implementation is much faster than ddm_rref for sparse matrices. In fact

1081 at the time of writing it is also (slightly) faster than the dense

1082 implementation even if the input is a fully dense matrix so it seems to be

1083 faster in all cases.

1084

1085 The elements of the matrix should support exact division with ``/``. For

1086 example elements of any domain that is a field (e.g. ``QQ``) should be

1087 fine. No attempt is made to handle inexact arithmetic.

1088

1089 """

1090 #

1091 # Any zeros in the matrix are not stored at all so an element is zero if

1092 # its row dict has no index at that key. A row is entirely zero if its

1093 # row index is not in the outer dict. Since rref reorders the rows and

1094 # removes zero rows we can completely discard the row indices. The first

1095 # step then copies the row dicts into a list sorted by the index of the

1096 # first nonzero column in each row.

1097 #

1098 # The algorithm then processes each row Ai one at a time. Previously seen

1099 # rows are used to cancel their pivot columns from Ai. Then a pivot from

1100 # Ai is chosen and is cancelled from all previously seen rows. At this

1101 # point Ai joins the previously seen rows. Once all rows are seen all

1102 # elimination has occurred and the rows are sorted by pivot column index.

1103 #

1104 # The previously seen rows are stored in two separate groups. The reduced

1105 # group consists of all rows that have been reduced to a single nonzero

1106 # element (the pivot). There is no need to attempt any further reduction

1107 # with these. Rows that still have other nonzeros need to be considered

1108 # when Ai is cancelled from the previously seen rows.

1109 #

1110 # A dict nonzerocolumns is used to map from a column index to a set of

1111 # previously seen rows that still have a nonzero element in that column.

1112 # This means that we can cancel the pivot from Ai into the previously seen

1113 # rows without needing to loop over each row that might have a zero in

1114 # that column.

1115 #

1116

1117 # Row dicts sorted by index of first nonzero column

1118 # (Maybe sorting is not needed/useful.)

1119 Arows = sorted((Ai.copy() for Ai in A.values()), key=min)

1120

1121 # Each processed row has an associated pivot column.

1122 # pivot_row_map maps from the pivot column index to the row dict.

1123 # This means that we can represent a set of rows purely as a set of their

1124 # pivot indices.

1125 pivot_row_map = {}

1126

1127 # Set of pivot indices for rows that are fully reduced to a single nonzero.

1128 reduced_pivots = set()

1129

1130 # Set of pivot indices for rows not fully reduced

1131 nonreduced_pivots = set()

1132

1133 # Map from column index to a set of pivot indices representing the rows

1134 # that have a nonzero at that column.

1135 nonzero_columns = defaultdict(set)

1136

1137 while Arows:

1138 # Select pivot element and row

1139 Ai = Arows.pop()

1140

1141 # Nonzero columns from fully reduced pivot rows can be removed

1142 Ai = {j: Aij for j, Aij in Ai.items() if j not in reduced_pivots}

1143

1144 # Others require full row cancellation

1145 for j in nonreduced_pivots & set(Ai):

1146 Aj = pivot_row_map[j]

1147 Aij = Ai[j]

1148 Ainz = set(Ai)

1149 Ajnz = set(Aj)

1150 for k in Ajnz - Ainz:

1151 Ai[k] = - Aij * Aj[k]

1152 Ai.pop(j)

1153 Ainz.remove(j)

1154 for k in Ajnz & Ainz:

1155 Aik = Ai[k] - Aij * Aj[k]

1156 if Aik:

1157 Ai[k] = Aik

1158 else:

1159 Ai.pop(k)

1160

1161 # We have now cancelled previously seen pivots from Ai.

1162 # If it is zero then discard it.

1163 if not Ai:

1164 continue

1165

1166 # Choose a pivot from Ai:

1167 j = min(Ai)

1168 Aij = Ai[j]

1169 pivot_row_map[j] = Ai

1170 Ainz = set(Ai)

1171

1172 # Normalise the pivot row to make the pivot 1.

1173 #

1174 # This approach is slow for some domains. Cross cancellation might be

1175 # better for e.g. QQ(x) with division delayed to the final steps.

1176 Aijinv = Aij**-1

1177 for l in Ai:

1178 Ai[l] *= Aijinv

1179

1180 # Use Aij to cancel column j from all previously seen rows

1181 for k in nonzero_columns.pop(j, ()):

1182 Ak = pivot_row_map[k]

1183 Akj = Ak[j]

1184 Aknz = set(Ak)

1185 for l in Ainz - Aknz:

1186 Ak[l] = - Akj * Ai[l]

1187 nonzero_columns[l].add(k)

1188 Ak.pop(j)

1189 Aknz.remove(j)

1190 for l in Ainz & Aknz:

1191 Akl = Ak[l] - Akj * Ai[l]

1192 if Akl:

1193 Ak[l] = Akl

1194 else:

1195 # Drop nonzero elements

1196 Ak.pop(l)

1197 if l != j:

1198 nonzero_columns[l].remove(k)

1199 if len(Ak) == 1:

1200 reduced_pivots.add(k)

1201 nonreduced_pivots.remove(k)

1202

1203 if len(Ai) == 1:

1204 reduced_pivots.add(j)

1205 else:

1206 nonreduced_pivots.add(j)

1207 for l in Ai:

1208 if l != j:

1209 nonzero_columns[l].add(j)

1210

1211 # All done!

1212 pivots = sorted(reduced_pivots | nonreduced_pivots)

1213 pivot2row = {p: n for n, p in enumerate(pivots)}

1214 nonzero_columns = {c: {pivot2row[p] for p in s} for c, s in nonzero_columns.items()}

1215 rows = [pivot_row_map[i] for i in pivots]

1216 rref = dict(enumerate(rows))

1217 return rref, pivots, nonzero_columns

1218

1219

1220def sdm_nullspace_from_rref(A, one, ncols, pivots, nonzero_cols):

1221 """Get nullspace from A which is in RREF"""

1222 nonpivots = sorted(set(range(ncols)) - set(pivots))

1223

1224 K = []

1225 for j in nonpivots:

1226 Kj = {j:one}

1227 for i in nonzero_cols.get(j, ()):

1228 Kj[pivots[i]] = -A[i][j]

1229 K.append(Kj)

1230

1231 return K, nonpivots

1232

1233

1234def sdm_particular_from_rref(A, ncols, pivots):

1235 """Get a particular solution from A which is in RREF"""

1236 P = {}

1237 for i, j in enumerate(pivots):

1238 Ain = A[i].get(ncols-1, None)

1239 if Ain is not None:

1240 P[j] = Ain / A[i][j]

1241 return P