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

1"""

3Module for the DDM class.

5The DDM class is an internal representation used by DomainMatrix. The letters

6DDM stand for Dense Domain Matrix. A DDM instance represents a matrix using

7elements from a polynomial Domain (e.g. ZZ, QQ, ...) in a dense-matrix

8representation.

10Basic usage:

12 >>> from sympy import ZZ, QQ

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

14 >>> A = DDM([[ZZ(0), ZZ(1)], [ZZ(-1), ZZ(0)]], (2, 2), ZZ)

15 >>> A.shape

16 (2, 2)

17 >>> A

18 [[0, 1], [-1, 0]]

19 >>> type(A)

20 <class 'sympy.polys.matrices.ddm.DDM'>

21 >>> A @ A

22 [[-1, 0], [0, -1]]

24The ddm_* functions are designed to operate on DDM as well as on an ordinary

25list of lists:

27 >>> from sympy.polys.matrices.dense import ddm_idet

28 >>> ddm_idet(A, QQ)

29 1

30 >>> ddm_idet([[0, 1], [-1, 0]], QQ)

31 1

32 >>> A

33 [[-1, 0], [0, -1]]

35Note that ddm_idet modifies the input matrix in-place. It is recommended to

36use the DDM.det method as a friendlier interface to this instead which takes

37care of copying the matrix:

39 >>> B = DDM([[ZZ(0), ZZ(1)], [ZZ(-1), ZZ(0)]], (2, 2), ZZ)

40 >>> B.det()

41 1

43Normally DDM would not be used directly and is just part of the internal

44representation of DomainMatrix which adds further functionality including e.g.

45unifying domains.

47The dense format used by DDM is a list of lists of elements e.g. the 2x2

48identity matrix is like [[1, 0], [0, 1]]. The DDM class itself is a subclass

49of list and its list items are plain lists. Elements are accessed as e.g.

50ddm[i][j] where ddm[i] gives the ith row and ddm[i][j] gets the element in the

51jth column of that row. Subclassing list makes e.g. iteration and indexing

52very efficient. We do not override __getitem__ because it would lose that

53benefit.

55The core routines are implemented by the ddm_* functions defined in dense.py.

56Those functions are intended to be able to operate on a raw list-of-lists

57representation of matrices with most functions operating in-place. The DDM

58class takes care of copying etc and also stores a Domain object associated

59with its elements. This makes it possible to implement things like A + B with

60domain checking and also shape checking so that the list of lists

61representation is friendlier.

63"""

64from itertools import chain

66from .exceptions import DMBadInputError, DMShapeError, DMDomainError

68from .dense import (

69 ddm_transpose,

70 ddm_iadd,

71 ddm_isub,

72 ddm_ineg,

73 ddm_imul,

74 ddm_irmul,

75 ddm_imatmul,

76 ddm_irref,

77 ddm_idet,

78 ddm_iinv,

79 ddm_ilu_split,

80 ddm_ilu_solve,

81 ddm_berk,

82 )

84from sympy.polys.domains import QQ

85from .lll import ddm_lll, ddm_lll_transform

88class DDM(list):

89 """Dense matrix based on polys domain elements

91 This is a list subclass and is a wrapper for a list of lists that supports

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

93 """

95 fmt = 'dense'

97 def __init__(self, rowslist, shape, domain):

98 super().__init__(rowslist)

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

100 self.domain = domain

101

102 if not (len(self) == m and all(len(row) == n for row in self)):

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

104

105 def getitem(self, i, j):

106 return self[i][j]

107

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

109 self[i][j] = value

110

111 def extract_slice(self, slice1, slice2):

112 ddm = [row[slice2] for row in self[slice1]]

113 rows = len(ddm)

114 cols = len(ddm[0]) if ddm else len(range(self.shape[1])[slice2])

115 return DDM(ddm, (rows, cols), self.domain)

116

117 def extract(self, rows, cols):

118 ddm = []

119 for i in rows:

120 rowi = self[i]

121 ddm.append([rowi[j] for j in cols])

122 return DDM(ddm, (len(rows), len(cols)), self.domain)

123

124 def to_list(self):

125 return list(self)

126

127 def to_list_flat(self):

128 flat = []

129 for row in self:

130 flat.extend(row)

131 return flat

132

133 def flatiter(self):

134 return chain.from_iterable(self)

135

136 def flat(self):

137 items = []

138 for row in self:

139 items.extend(row)

140 return items

141

142 def to_dok(self):

143 return {(i, j): e for i, row in enumerate(self) for j, e in enumerate(row)}

144

145 def to_ddm(self):

146 return self

147

148 def to_sdm(self):

149 return SDM.from_list(self, self.shape, self.domain)

150

151 def convert_to(self, K):

152 Kold = self.domain

153 if K == Kold:

154 return self.copy()

155 rows = ([K.convert_from(e, Kold) for e in row] for row in self)

156 return DDM(rows, self.shape, K)

157

158 def __str__(self):

159 rowsstr = ['[%s]' % ', '.join(map(str, row)) for row in self]

160 return '[%s]' % ', '.join(rowsstr)

161

162 def __repr__(self):

163 cls = type(self).__name__

164 rows = list.__repr__(self)

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

166

167 def __eq__(self, other):

168 if not isinstance(other, DDM):

169 return False

170 return (super().__eq__(other) and self.domain == other.domain)

171

172 def __ne__(self, other):

173 return not self.__eq__(other)

174

175 @classmethod

176 def zeros(cls, shape, domain):

177 z = domain.zero

178 m, n = shape

179 rowslist = ([z] * n for _ in range(m))

180 return DDM(rowslist, shape, domain)

181

182 @classmethod

183 def ones(cls, shape, domain):

184 one = domain.one

185 m, n = shape

186 rowlist = ([one] * n for _ in range(m))

187 return DDM(rowlist, shape, domain)

188

189 @classmethod

190 def eye(cls, size, domain):

191 one = domain.one

192 ddm = cls.zeros((size, size), domain)

193 for i in range(size):

194 ddm[i][i] = one

195 return ddm

196

197 def copy(self):

198 copyrows = (row[:] for row in self)

199 return DDM(copyrows, self.shape, self.domain)

200

201 def transpose(self):

202 rows, cols = self.shape

203 if rows:

204 ddmT = ddm_transpose(self)

205 else:

206 ddmT = [[]] * cols

207 return DDM(ddmT, (cols, rows), self.domain)

208

209 def __add__(a, b):

210 if not isinstance(b, DDM):

211 return NotImplemented

212 return a.add(b)

213

214 def __sub__(a, b):

215 if not isinstance(b, DDM):

216 return NotImplemented

217 return a.sub(b)

218

219 def __neg__(a):

220 return a.neg()

221

222 def __mul__(a, b):

223 if b in a.domain:

224 return a.mul(b)

225 else:

226 return NotImplemented

227

228 def __rmul__(a, b):

229 if b in a.domain:

230 return a.mul(b)

231 else:

232 return NotImplemented

233

234 def __matmul__(a, b):

235 if isinstance(b, DDM):

236 return a.matmul(b)

237 else:

238 return NotImplemented

239

240 @classmethod

241 def _check(cls, a, op, b, ashape, bshape):

242 if a.domain != b.domain:

243 msg = "Domain mismatch: %s %s %s" % (a.domain, op, b.domain)

244 raise DMDomainError(msg)

245 if ashape != bshape:

246 msg = "Shape mismatch: %s %s %s" % (a.shape, op, b.shape)

247 raise DMShapeError(msg)

248

249 def add(a, b):

250 """a + b"""

251 a._check(a, '+', b, a.shape, b.shape)

252 c = a.copy()

253 ddm_iadd(c, b)

254 return c

255

256 def sub(a, b):

257 """a - b"""

258 a._check(a, '-', b, a.shape, b.shape)

259 c = a.copy()

260 ddm_isub(c, b)

261 return c

262

263 def neg(a):

264 """-a"""

265 b = a.copy()

266 ddm_ineg(b)

267 return b

268

269 def mul(a, b):

270 c = a.copy()

271 ddm_imul(c, b)

272 return c

273

274 def rmul(a, b):

275 c = a.copy()

276 ddm_irmul(c, b)

277 return c

278

279 def matmul(a, b):

280 """a @ b (matrix product)"""

281 m, o = a.shape

282 o2, n = b.shape

283 a._check(a, '*', b, o, o2)

284 c = a.zeros((m, n), a.domain)

285 ddm_imatmul(c, a, b)

286 return c

287

288 def mul_elementwise(a, b):

289 assert a.shape == b.shape

290 assert a.domain == b.domain

291 c = [[aij * bij for aij, bij in zip(ai, bi)] for ai, bi in zip(a, b)]

292 return DDM(c, a.shape, a.domain)

293

294 def hstack(A, *B):

295 """Horizontally stacks :py:class:`~.DDM` matrices.

296

297 Examples

298 ========

299

300 >>> from sympy import ZZ

301 >>> from sympy.polys.matrices.sdm import DDM

302

303 >>> A = DDM([[ZZ(1), ZZ(2)], [ZZ(3), ZZ(4)]], (2, 2), ZZ)

304 >>> B = DDM([[ZZ(5), ZZ(6)], [ZZ(7), ZZ(8)]], (2, 2), ZZ)

305 >>> A.hstack(B)

306 [[1, 2, 5, 6], [3, 4, 7, 8]]

307

308 >>> C = DDM([[ZZ(9), ZZ(10)], [ZZ(11), ZZ(12)]], (2, 2), ZZ)

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

310 [[1, 2, 5, 6, 9, 10], [3, 4, 7, 8, 11, 12]]

311 """

312 Anew = list(A.copy())

313 rows, cols = A.shape

314 domain = A.domain

315

316 for Bk in B:

317 Bkrows, Bkcols = Bk.shape

318 assert Bkrows == rows

319 assert Bk.domain == domain

320

321 cols += Bkcols

322

323 for i, Bki in enumerate(Bk):

324 Anew[i].extend(Bki)

325

326 return DDM(Anew, (rows, cols), A.domain)

327

328 def vstack(A, *B):

329 """Vertically stacks :py:class:`~.DDM` matrices.

330

331 Examples

332 ========

333

334 >>> from sympy import ZZ

335 >>> from sympy.polys.matrices.sdm import DDM

336

337 >>> A = DDM([[ZZ(1), ZZ(2)], [ZZ(3), ZZ(4)]], (2, 2), ZZ)

338 >>> B = DDM([[ZZ(5), ZZ(6)], [ZZ(7), ZZ(8)]], (2, 2), ZZ)

339 >>> A.vstack(B)

340 [[1, 2], [3, 4], [5, 6], [7, 8]]

341

342 >>> C = DDM([[ZZ(9), ZZ(10)], [ZZ(11), ZZ(12)]], (2, 2), ZZ)

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

344 [[1, 2], [3, 4], [5, 6], [7, 8], [9, 10], [11, 12]]

345 """

346 Anew = list(A.copy())

347 rows, cols = A.shape

348 domain = A.domain

349

350 for Bk in B:

351 Bkrows, Bkcols = Bk.shape

352 assert Bkcols == cols

353 assert Bk.domain == domain

354

355 rows += Bkrows

356

357 Anew.extend(Bk.copy())

358

359 return DDM(Anew, (rows, cols), A.domain)

360

361 def applyfunc(self, func, domain):

362 elements = (list(map(func, row)) for row in self)

363 return DDM(elements, self.shape, domain)

364

365 def scc(a):

366 """Strongly connected components of a square matrix *a*.

367

368 Examples

369 ========

370

371 >>> from sympy import ZZ

372 >>> from sympy.polys.matrices.sdm import DDM

373 >>> A = DDM([[ZZ(1), ZZ(0)], [ZZ(0), ZZ(1)]], (2, 2), ZZ)

374 >>> A.scc()

375 [[0], [1]]

376

377 See also

378 ========

379

380 sympy.polys.matrices.domainmatrix.DomainMatrix.scc

381

382 """

383 return a.to_sdm().scc()

384

385 def rref(a):

386 """Reduced-row echelon form of a and list of pivots"""

387 b = a.copy()

388 K = a.domain

389 partial_pivot = K.is_RealField or K.is_ComplexField

390 pivots = ddm_irref(b, _partial_pivot=partial_pivot)

391 return b, pivots

392

393 def nullspace(a):

394 rref, pivots = a.rref()

395 rows, cols = a.shape

396 domain = a.domain

397

398 basis = []

399 nonpivots = []

400 for i in range(cols):

401 if i in pivots:

402 continue

403 nonpivots.append(i)

404 vec = [domain.one if i == j else domain.zero for j in range(cols)]

405 for ii, jj in enumerate(pivots):

406 vec[jj] -= rref[ii][i]

407 basis.append(vec)

408

409 return DDM(basis, (len(basis), cols), domain), nonpivots

410

411 def particular(a):

412 return a.to_sdm().particular().to_ddm()

413

414 def det(a):

415 """Determinant of a"""

416 m, n = a.shape

417 if m != n:

418 raise DMShapeError("Determinant of non-square matrix")

419 b = a.copy()

420 K = b.domain

421 deta = ddm_idet(b, K)

422 return deta

423

424 def inv(a):

425 """Inverse of a"""

426 m, n = a.shape

427 if m != n:

428 raise DMShapeError("Determinant of non-square matrix")

429 ainv = a.copy()

430 K = a.domain

431 ddm_iinv(ainv, a, K)

432 return ainv

433

434 def lu(a):

435 """L, U decomposition of a"""

436 m, n = a.shape

437 K = a.domain

438

439 U = a.copy()

440 L = a.eye(m, K)

441 swaps = ddm_ilu_split(L, U, K)

442

443 return L, U, swaps

444

445 def lu_solve(a, b):

446 """x where a*x = b"""

447 m, n = a.shape

448 m2, o = b.shape

449 a._check(a, 'lu_solve', b, m, m2)

450

451 L, U, swaps = a.lu()

452 x = a.zeros((n, o), a.domain)

453 ddm_ilu_solve(x, L, U, swaps, b)

454 return x

455

456 def charpoly(a):

457 """Coefficients of characteristic polynomial of a"""

458 K = a.domain

459 m, n = a.shape

460 if m != n:

461 raise DMShapeError("Charpoly of non-square matrix")

462 vec = ddm_berk(a, K)

463 coeffs = [vec[i][0] for i in range(n+1)]

464 return coeffs

465

466 def is_zero_matrix(self):

467 """

468 Says whether this matrix has all zero entries.

469 """

470 zero = self.domain.zero

471 return all(Mij == zero for Mij in self.flatiter())

472

473 def is_upper(self):

474 """

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

476 even if the matrix is not square.

477 """

478 zero = self.domain.zero

479 return all(Mij == zero for i, Mi in enumerate(self) for Mij in Mi[:i])

480

481 def is_lower(self):

482 """

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

484 even if the matrix is not square.

485 """

486 zero = self.domain.zero

487 return all(Mij == zero for i, Mi in enumerate(self) for Mij in Mi[i+1:])

488

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

490 return ddm_lll(A, delta=delta)

491

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

493 return ddm_lll_transform(A, delta=delta)

494

495

496from .sdm import SDM