Coverage for /usr/lib/python3/dist-packages/sympy/matrices/expressions/matexpr.py: 39%

1from __future__ import annotations 

2from functools import wraps 

3 

4from sympy.core import S, Integer, Basic, Mul, Add 

5from sympy.core.assumptions import check_assumptions 

6from sympy.core.decorators import call_highest_priority 

7from sympy.core.expr import Expr, ExprBuilder 

8from sympy.core.logic import FuzzyBool 

9from sympy.core.symbol import Str, Dummy, symbols, Symbol 

10from sympy.core.sympify import SympifyError, _sympify 

11from sympy.external.gmpy import SYMPY_INTS 

12from sympy.functions import conjugate, adjoint 

13from sympy.functions.special.tensor_functions import KroneckerDelta 

14from sympy.matrices.common import NonSquareMatrixError 

15from sympy.matrices.matrices import MatrixKind, MatrixBase 

16from sympy.multipledispatch import dispatch 

17from sympy.utilities.misc import filldedent 

18 

19 

20def _sympifyit(arg, retval=None): 

21 # This version of _sympifyit sympifies MutableMatrix objects 

22 def deco(func): 

23 @wraps(func) 

24 def __sympifyit_wrapper(a, b): 

25 try: 

26 b = _sympify(b) 

27 return func(a, b) 

28 except SympifyError: 

29 return retval 

30 

31 return __sympifyit_wrapper 

32 

33 return deco 

34 

35 

36class MatrixExpr(Expr): 

37 """Superclass for Matrix Expressions 

38 

39 MatrixExprs represent abstract matrices, linear transformations represented 

40 within a particular basis. 

41 

42 Examples 

43 ======== 

44 

45 >>> from sympy import MatrixSymbol 

46 >>> A = MatrixSymbol('A', 3, 3) 

47 >>> y = MatrixSymbol('y', 3, 1) 

48 >>> x = (A.T*A).I * A * y 

49 

50 See Also 

51 ======== 

52 

53 MatrixSymbol, MatAdd, MatMul, Transpose, Inverse 

54 """ 

55 __slots__: tuple[str, ...] = () 

56 

57 # Should not be considered iterable by the 

58 # sympy.utilities.iterables.iterable function. Subclass that actually are 

59 # iterable (i.e., explicit matrices) should set this to True. 

60 _iterable = False 

61 

62 _op_priority = 11.0 

63 

64 is_Matrix: bool = True 

65 is_MatrixExpr: bool = True 

66 is_Identity: FuzzyBool = None 

67 is_Inverse = False 

68 is_Transpose = False 

69 is_ZeroMatrix = False 

70 is_MatAdd = False 

71 is_MatMul = False 

72 

73 is_commutative = False 

74 is_number = False 

75 is_symbol = False 

76 is_scalar = False 

77 

78 kind: MatrixKind = MatrixKind() 

79 

80 def __new__(cls, *args, **kwargs): 

81 args = map(_sympify, args) 

82 return Basic.__new__(cls, *args, **kwargs) 

83 

84 # The following is adapted from the core Expr object 

85 

86 @property 

87 def shape(self) -> tuple[Expr, Expr]: 

88 raise NotImplementedError 

89 

90 @property 

91 def _add_handler(self): 

92 return MatAdd 

93 

94 @property 

95 def _mul_handler(self): 

96 return MatMul 

97 

98 def __neg__(self): 

99 return MatMul(S.NegativeOne, self).doit() 

100 

101 def __abs__(self): 

102 raise NotImplementedError 

103 

104 @_sympifyit('other', NotImplemented) 

105 @call_highest_priority('__radd__') 

106 def __add__(self, other): 

107 return MatAdd(self, other).doit() 

108 

109 @_sympifyit('other', NotImplemented) 

110 @call_highest_priority('__add__') 

111 def __radd__(self, other): 

112 return MatAdd(other, self).doit() 

113 

114 @_sympifyit('other', NotImplemented) 

115 @call_highest_priority('__rsub__') 

116 def __sub__(self, other): 

117 return MatAdd(self, -other).doit() 

118 

119 @_sympifyit('other', NotImplemented) 

120 @call_highest_priority('__sub__') 

121 def __rsub__(self, other): 

122 return MatAdd(other, -self).doit() 

123 

124 @_sympifyit('other', NotImplemented) 

125 @call_highest_priority('__rmul__') 

126 def __mul__(self, other): 

127 return MatMul(self, other).doit() 

128 

129 @_sympifyit('other', NotImplemented) 

130 @call_highest_priority('__rmul__') 

131 def __matmul__(self, other): 

132 return MatMul(self, other).doit() 

133 

134 @_sympifyit('other', NotImplemented) 

135 @call_highest_priority('__mul__') 

136 def __rmul__(self, other): 

137 return MatMul(other, self).doit() 

138 

139 @_sympifyit('other', NotImplemented) 

140 @call_highest_priority('__mul__') 

141 def __rmatmul__(self, other): 

142 return MatMul(other, self).doit() 

143 

144 @_sympifyit('other', NotImplemented) 

145 @call_highest_priority('__rpow__') 

146 def __pow__(self, other): 

147 return MatPow(self, other).doit() 

148 

149 @_sympifyit('other', NotImplemented) 

150 @call_highest_priority('__pow__') 

151 def __rpow__(self, other): 

152 raise NotImplementedError("Matrix Power not defined") 

153 

154 @_sympifyit('other', NotImplemented) 

155 @call_highest_priority('__rtruediv__') 

156 def __truediv__(self, other): 

157 return self * other**S.NegativeOne 

158 

159 @_sympifyit('other', NotImplemented) 

160 @call_highest_priority('__truediv__') 

161 def __rtruediv__(self, other): 

162 raise NotImplementedError() 

163 #return MatMul(other, Pow(self, S.NegativeOne)) 

164 

165 @property 

166 def rows(self): 

167 return self.shape[0] 

168 

169 @property 

170 def cols(self): 

171 return self.shape[1] 

172 

173 @property 

174 def is_square(self) -> bool | None: 

175 rows, cols = self.shape 

176 if isinstance(rows, Integer) and isinstance(cols, Integer): 

177 return rows == cols 

178 if rows == cols: 

179 return True 

180 return None 

181 

182 def _eval_conjugate(self): 

183 from sympy.matrices.expressions.adjoint import Adjoint 

184 return Adjoint(Transpose(self)) 

185 

186 def as_real_imag(self, deep=True, **hints): 

187 return self._eval_as_real_imag() 

188 

189 def _eval_as_real_imag(self): 

190 real = S.Half * (self + self._eval_conjugate()) 

191 im = (self - self._eval_conjugate())/(2*S.ImaginaryUnit) 

192 return (real, im) 

193 

194 def _eval_inverse(self): 

195 return Inverse(self) 

196 

197 def _eval_determinant(self): 

198 return Determinant(self) 

199 

200 def _eval_transpose(self): 

201 return Transpose(self) 

202 

203 def _eval_power(self, exp): 

204 """ 

205 Override this in sub-classes to implement simplification of powers. The cases where the exponent 

206 is -1, 0, 1 are already covered in MatPow.doit(), so implementations can exclude these cases. 

207 """ 

208 return MatPow(self, exp) 

209 

210 def _eval_simplify(self, **kwargs): 

211 if self.is_Atom: 

212 return self 

213 else: 

214 from sympy.simplify import simplify 

215 return self.func(*[simplify(x, **kwargs) for x in self.args]) 

216 

217 def _eval_adjoint(self): 

218 from sympy.matrices.expressions.adjoint import Adjoint 

219 return Adjoint(self) 

220 

221 def _eval_derivative_n_times(self, x, n): 

222 return Basic._eval_derivative_n_times(self, x, n) 

223 

224 def _eval_derivative(self, x): 

225 # `x` is a scalar: 

226 if self.has(x): 

227 # See if there are other methods using it: 

228 return super()._eval_derivative(x) 

229 else: 

230 return ZeroMatrix(*self.shape) 

231 

232 @classmethod 

233 def _check_dim(cls, dim): 

234 """Helper function to check invalid matrix dimensions""" 

235 ok = check_assumptions(dim, integer=True, nonnegative=True) 

236 if ok is False: 

237 raise ValueError( 

238 "The dimension specification {} should be " 

239 "a nonnegative integer.".format(dim)) 

240 

241 

242 def _entry(self, i, j, **kwargs): 

243 raise NotImplementedError( 

244 "Indexing not implemented for %s" % self.__class__.__name__) 

245 

246 def adjoint(self): 

247 return adjoint(self) 

248 

249 def as_coeff_Mul(self, rational=False): 

250 """Efficiently extract the coefficient of a product.""" 

251 return S.One, self 

252 

253 def conjugate(self): 

254 return conjugate(self) 

255 

256 def transpose(self): 

257 from sympy.matrices.expressions.transpose import transpose 

258 return transpose(self) 

259 

260 @property 

261 def T(self): 

262 '''Matrix transposition''' 

263 return self.transpose() 

264 

265 def inverse(self): 

266 if self.is_square is False: 

267 raise NonSquareMatrixError('Inverse of non-square matrix') 

268 return self._eval_inverse() 

269 

270 def inv(self): 

271 return self.inverse() 

272 

273 def det(self): 

274 from sympy.matrices.expressions.determinant import det 

275 return det(self) 

276 

277 @property 

278 def I(self): 

279 return self.inverse() 

280 

281 def valid_index(self, i, j): 

282 def is_valid(idx): 

283 return isinstance(idx, (int, Integer, Symbol, Expr)) 

284 return (is_valid(i) and is_valid(j) and 

285 (self.rows is None or 

286 (i >= -self.rows) != False and (i < self.rows) != False) and 

287 (j >= -self.cols) != False and (j < self.cols) != False) 

288 

289 def __getitem__(self, key): 

290 if not isinstance(key, tuple) and isinstance(key, slice): 

291 from sympy.matrices.expressions.slice import MatrixSlice 

292 return MatrixSlice(self, key, (0, None, 1)) 

293 if isinstance(key, tuple) and len(key) == 2: 

294 i, j = key 

295 if isinstance(i, slice) or isinstance(j, slice): 

296 from sympy.matrices.expressions.slice import MatrixSlice 

297 return MatrixSlice(self, i, j) 

298 i, j = _sympify(i), _sympify(j) 

299 if self.valid_index(i, j) != False: 

300 return self._entry(i, j) 

301 else: 

302 raise IndexError("Invalid indices (%s, %s)" % (i, j)) 

303 elif isinstance(key, (SYMPY_INTS, Integer)): 

304 # row-wise decomposition of matrix 

305 rows, cols = self.shape 

306 # allow single indexing if number of columns is known 

307 if not isinstance(cols, Integer): 

308 raise IndexError(filldedent(''' 

309 Single indexing is only supported when the number 

310 of columns is known.''')) 

311 key = _sympify(key) 

312 i = key // cols 

313 j = key % cols 

314 if self.valid_index(i, j) != False: 

315 return self._entry(i, j) 

316 else: 

317 raise IndexError("Invalid index %s" % key) 

318 elif isinstance(key, (Symbol, Expr)): 

319 raise IndexError(filldedent(''' 

320 Only integers may be used when addressing the matrix 

321 with a single index.''')) 

322 raise IndexError("Invalid index, wanted %s[i,j]" % self) 

323 

324 def _is_shape_symbolic(self) -> bool: 

325 return (not isinstance(self.rows, (SYMPY_INTS, Integer)) 

326 or not isinstance(self.cols, (SYMPY_INTS, Integer))) 

327 

328 def as_explicit(self): 

329 """ 

330 Returns a dense Matrix with elements represented explicitly 

331 

332 Returns an object of type ImmutableDenseMatrix. 

333 

334 Examples 

335 ======== 

336 

337 >>> from sympy import Identity 

338 >>> I = Identity(3) 

339 >>> I 

340 I 

341 >>> I.as_explicit() 

342 Matrix([ 

343 [1, 0, 0], 

344 [0, 1, 0], 

345 [0, 0, 1]]) 

346 

347 See Also 

348 ======== 

349 as_mutable: returns mutable Matrix type 

350 

351 """ 

352 if self._is_shape_symbolic(): 

353 raise ValueError( 

354 'Matrix with symbolic shape ' 

355 'cannot be represented explicitly.') 

356 from sympy.matrices.immutable import ImmutableDenseMatrix 

357 return ImmutableDenseMatrix([[self[i, j] 

358 for j in range(self.cols)] 

359 for i in range(self.rows)]) 

360 

361 def as_mutable(self): 

362 """ 

363 Returns a dense, mutable matrix with elements represented explicitly 

364 

365 Examples 

366 ======== 

367 

368 >>> from sympy import Identity 

369 >>> I = Identity(3) 

370 >>> I 

371 I 

372 >>> I.shape 

373 (3, 3) 

374 >>> I.as_mutable() 

375 Matrix([ 

376 [1, 0, 0], 

377 [0, 1, 0], 

378 [0, 0, 1]]) 

379 

380 See Also 

381 ======== 

382 as_explicit: returns ImmutableDenseMatrix 

383 """ 

384 return self.as_explicit().as_mutable() 

385 

386 def __array__(self): 

387 from numpy import empty 

388 a = empty(self.shape, dtype=object) 

389 for i in range(self.rows): 

390 for j in range(self.cols): 

391 a[i, j] = self[i, j] 

392 return a 

393 

394 def equals(self, other): 

395 """ 

396 Test elementwise equality between matrices, potentially of different 

397 types 

398 

399 >>> from sympy import Identity, eye 

400 >>> Identity(3).equals(eye(3)) 

401 True 

402 """ 

403 return self.as_explicit().equals(other) 

404 

405 def canonicalize(self): 

406 return self 

407 

408 def as_coeff_mmul(self): 

409 return S.One, MatMul(self) 

410 

411 @staticmethod 

412 def from_index_summation(expr, first_index=None, last_index=None, dimensions=None): 

413 r""" 

414 Parse expression of matrices with explicitly summed indices into a 

415 matrix expression without indices, if possible. 

416 

417 This transformation expressed in mathematical notation: 

418 

419 `\sum_{j=0}^{N-1} A_{i,j} B_{j,k} \Longrightarrow \mathbf{A}\cdot \mathbf{B}` 

420 

421 Optional parameter ``first_index``: specify which free index to use as 

422 the index starting the expression. 

423 

424 Examples 

425 ======== 

426 

427 >>> from sympy import MatrixSymbol, MatrixExpr, Sum 

428 >>> from sympy.abc import i, j, k, l, N 

429 >>> A = MatrixSymbol("A", N, N) 

430 >>> B = MatrixSymbol("B", N, N) 

431 >>> expr = Sum(A[i, j]*B[j, k], (j, 0, N-1)) 

432 >>> MatrixExpr.from_index_summation(expr) 

433 A*B 

434 

435 Transposition is detected: 

436 

437 >>> expr = Sum(A[j, i]*B[j, k], (j, 0, N-1)) 

438 >>> MatrixExpr.from_index_summation(expr) 

439 A.T*B 

440 

441 Detect the trace: 

442 

443 >>> expr = Sum(A[i, i], (i, 0, N-1)) 

444 >>> MatrixExpr.from_index_summation(expr) 

445 Trace(A) 

446 

447 More complicated expressions: 

448 

449 >>> expr = Sum(A[i, j]*B[k, j]*A[l, k], (j, 0, N-1), (k, 0, N-1)) 

450 >>> MatrixExpr.from_index_summation(expr) 

451 A*B.T*A.T 

452 """ 

453 from sympy.tensor.array.expressions.from_indexed_to_array import convert_indexed_to_array 

454 from sympy.tensor.array.expressions.from_array_to_matrix import convert_array_to_matrix 

455 first_indices = [] 

456 if first_index is not None: 

457 first_indices.append(first_index) 

458 if last_index is not None: 

459 first_indices.append(last_index) 

460 arr = convert_indexed_to_array(expr, first_indices=first_indices) 

461 return convert_array_to_matrix(arr) 

462 

463 def applyfunc(self, func): 

464 from .applyfunc import ElementwiseApplyFunction 

465 return ElementwiseApplyFunction(func, self) 

466 

467 

468@dispatch(MatrixExpr, Expr) 

469def _eval_is_eq(lhs, rhs): # noqa:F811 

470 return False 

471 

472@dispatch(MatrixExpr, MatrixExpr) # type: ignore 

473def _eval_is_eq(lhs, rhs): # noqa:F811 

474 if lhs.shape != rhs.shape: 

475 return False 

476 if (lhs - rhs).is_ZeroMatrix: 

477 return True 

478 

479def get_postprocessor(cls): 

480 def _postprocessor(expr): 

481 # To avoid circular imports, we can't have MatMul/MatAdd on the top level 

482 mat_class = {Mul: MatMul, Add: MatAdd}[cls] 

483 nonmatrices = [] 

484 matrices = [] 

485 for term in expr.args: 

486 if isinstance(term, MatrixExpr): 

487 matrices.append(term) 

488 else: 

489 nonmatrices.append(term) 

490 

491 if not matrices: 

492 return cls._from_args(nonmatrices) 

493 

494 if nonmatrices: 

495 if cls == Mul: 

496 for i in range(len(matrices)): 

497 if not matrices[i].is_MatrixExpr: 

498 # If one of the matrices explicit, absorb the scalar into it 

499 # (doit will combine all explicit matrices into one, so it 

500 # doesn't matter which) 

501 matrices[i] = matrices[i].__mul__(cls._from_args(nonmatrices)) 

502 nonmatrices = [] 

503 break 

504 

505 else: 

506 # Maintain the ability to create Add(scalar, matrix) without 

507 # raising an exception. That way different algorithms can 

508 # replace matrix expressions with non-commutative symbols to 

509 # manipulate them like non-commutative scalars. 

510 return cls._from_args(nonmatrices + [mat_class(*matrices).doit(deep=False)]) 

511 

512 if mat_class == MatAdd: 

513 return mat_class(*matrices).doit(deep=False) 

514 return mat_class(cls._from_args(nonmatrices), *matrices).doit(deep=False) 

515 return _postprocessor 

516 

517 

518Basic._constructor_postprocessor_mapping[MatrixExpr] = { 

519 "Mul": [get_postprocessor(Mul)], 

520 "Add": [get_postprocessor(Add)], 

521} 

522 

523 

524def _matrix_derivative(expr, x, old_algorithm=False): 

525 

526 if isinstance(expr, MatrixBase) or isinstance(x, MatrixBase): 

527 # Do not use array expressions for explicit matrices: 

528 old_algorithm = True 

529 

530 if old_algorithm: 

531 return _matrix_derivative_old_algorithm(expr, x) 

532 

533 from sympy.tensor.array.expressions.from_matrix_to_array import convert_matrix_to_array 

534 from sympy.tensor.array.expressions.arrayexpr_derivatives import array_derive 

535 from sympy.tensor.array.expressions.from_array_to_matrix import convert_array_to_matrix 

536 

537 array_expr = convert_matrix_to_array(expr) 

538 diff_array_expr = array_derive(array_expr, x) 

539 diff_matrix_expr = convert_array_to_matrix(diff_array_expr) 

540 return diff_matrix_expr 

541 

542 

543def _matrix_derivative_old_algorithm(expr, x): 

544 from sympy.tensor.array.array_derivatives import ArrayDerivative 

545 lines = expr._eval_derivative_matrix_lines(x) 

546 

547 parts = [i.build() for i in lines] 

548 

549 from sympy.tensor.array.expressions.from_array_to_matrix import convert_array_to_matrix 

550 

551 parts = [[convert_array_to_matrix(j) for j in i] for i in parts] 

552 

553 def _get_shape(elem): 

554 if isinstance(elem, MatrixExpr): 

555 return elem.shape 

556 return 1, 1 

557 

558 def get_rank(parts): 

559 return sum([j not in (1, None) for i in parts for j in _get_shape(i)]) 

560 

561 ranks = [get_rank(i) for i in parts] 

562 rank = ranks[0] 

563 

564 def contract_one_dims(parts): 

565 if len(parts) == 1: 

566 return parts[0] 

567 else: 

568 p1, p2 = parts[:2] 

569 if p2.is_Matrix: 

570 p2 = p2.T 

571 if p1 == Identity(1): 

572 pbase = p2 

573 elif p2 == Identity(1): 

574 pbase = p1 

575 else: 

576 pbase = p1*p2 

577 if len(parts) == 2: 

578 return pbase 

579 else: # len(parts) > 2 

580 if pbase.is_Matrix: 

581 raise ValueError("") 

582 return pbase*Mul.fromiter(parts[2:]) 

583 

584 if rank <= 2: 

585 return Add.fromiter([contract_one_dims(i) for i in parts]) 

586 

587 return ArrayDerivative(expr, x) 

588 

589 

590class MatrixElement(Expr): 

591 parent = property(lambda self: self.args[0]) 

592 i = property(lambda self: self.args[1]) 

593 j = property(lambda self: self.args[2]) 

594 _diff_wrt = True 

595 is_symbol = True 

596 is_commutative = True 

597 

598 def __new__(cls, name, n, m): 

599 n, m = map(_sympify, (n, m)) 

600 from sympy.matrices.matrices import MatrixBase 

601 if isinstance(name, str): 

602 name = Symbol(name) 

603 else: 

604 if isinstance(name, MatrixBase): 

605 if n.is_Integer and m.is_Integer: 

606 return name[n, m] 

607 name = _sympify(name) # change mutable into immutable 

608 else: 

609 name = _sympify(name) 

610 if not isinstance(name.kind, MatrixKind): 

611 raise TypeError("First argument of MatrixElement should be a matrix") 

612 if not getattr(name, 'valid_index', lambda n, m: True)(n, m): 

613 raise IndexError('indices out of range') 

614 obj = Expr.__new__(cls, name, n, m) 

615 return obj 

616 

617 @property 

618 def symbol(self): 

619 return self.args[0] 

620 

621 def doit(self, **hints): 

622 deep = hints.get('deep', True) 

623 if deep: 

624 args = [arg.doit(**hints) for arg in self.args] 

625 else: 

626 args = self.args 

627 return args[0][args[1], args[2]] 

628 

629 @property 

630 def indices(self): 

631 return self.args[1:] 

632 

633 def _eval_derivative(self, v): 

634 

635 if not isinstance(v, MatrixElement): 

636 from sympy.matrices.matrices import MatrixBase 

637 if isinstance(self.parent, MatrixBase): 

638 return self.parent.diff(v)[self.i, self.j] 

639 return S.Zero 

640 

641 M = self.args[0] 

642 

643 m, n = self.parent.shape 

644 

645 if M == v.args[0]: 

646 return KroneckerDelta(self.args[1], v.args[1], (0, m-1)) * \ 

647 KroneckerDelta(self.args[2], v.args[2], (0, n-1)) 

648 

649 if isinstance(M, Inverse): 

650 from sympy.concrete.summations import Sum 

651 i, j = self.args[1:] 

652 i1, i2 = symbols("z1, z2", cls=Dummy) 

653 Y = M.args[0] 

654 r1, r2 = Y.shape 

655 return -Sum(M[i, i1]*Y[i1, i2].diff(v)*M[i2, j], (i1, 0, r1-1), (i2, 0, r2-1)) 

656 

657 if self.has(v.args[0]): 

658 return None 

659 

660 return S.Zero 

661 

662 

663class MatrixSymbol(MatrixExpr): 

664 """Symbolic representation of a Matrix object 

665 

666 Creates a SymPy Symbol to represent a Matrix. This matrix has a shape and 

667 can be included in Matrix Expressions 

668 

669 Examples 

670 ======== 

671 

672 >>> from sympy import MatrixSymbol, Identity 

673 >>> A = MatrixSymbol('A', 3, 4) # A 3 by 4 Matrix 

674 >>> B = MatrixSymbol('B', 4, 3) # A 4 by 3 Matrix 

675 >>> A.shape 

676 (3, 4) 

677 >>> 2*A*B + Identity(3) 

678 I + 2*A*B 

679 """ 

680 is_commutative = False 

681 is_symbol = True 

682 _diff_wrt = True 

683 

684 def __new__(cls, name, n, m): 

685 n, m = _sympify(n), _sympify(m) 

686 

687 cls._check_dim(m) 

688 cls._check_dim(n) 

689 

690 if isinstance(name, str): 

691 name = Str(name) 

692 obj = Basic.__new__(cls, name, n, m) 

693 return obj 

694 

695 @property 

696 def shape(self): 

697 return self.args[1], self.args[2] 

698 

699 @property 

700 def name(self): 

701 return self.args[0].name 

702 

703 def _entry(self, i, j, **kwargs): 

704 return MatrixElement(self, i, j) 

705 

706 @property 

707 def free_symbols(self): 

708 return {self} 

709 

710 def _eval_simplify(self, **kwargs): 

711 return self 

712 

713 def _eval_derivative(self, x): 

714 # x is a scalar: 

715 return ZeroMatrix(self.shape[0], self.shape[1]) 

716 

717 def _eval_derivative_matrix_lines(self, x): 

718 if self != x: 

719 first = ZeroMatrix(x.shape[0], self.shape[0]) if self.shape[0] != 1 else S.Zero 

720 second = ZeroMatrix(x.shape[1], self.shape[1]) if self.shape[1] != 1 else S.Zero 

721 return [_LeftRightArgs( 

722 [first, second], 

723 )] 

724 else: 

725 first = Identity(self.shape[0]) if self.shape[0] != 1 else S.One 

726 second = Identity(self.shape[1]) if self.shape[1] != 1 else S.One 

727 return [_LeftRightArgs( 

728 [first, second], 

729 )] 

730 

731 

732def matrix_symbols(expr): 

733 return [sym for sym in expr.free_symbols if sym.is_Matrix] 

734 

735 

736class _LeftRightArgs: 

737 r""" 

738 Helper class to compute matrix derivatives. 

739 

740 The logic: when an expression is derived by a matrix `X_{mn}`, two lines of 

741 matrix multiplications are created: the one contracted to `m` (first line), 

742 and the one contracted to `n` (second line). 

743 

744 Transposition flips the side by which new matrices are connected to the 

745 lines. 

746 

747 The trace connects the end of the two lines. 

748 """ 

749 

750 def __init__(self, lines, higher=S.One): 

751 self._lines = list(lines) 

752 self._first_pointer_parent = self._lines 

753 self._first_pointer_index = 0 

754 self._first_line_index = 0 

755 self._second_pointer_parent = self._lines 

756 self._second_pointer_index = 1 

757 self._second_line_index = 1 

758 self.higher = higher 

759 

760 @property 

761 def first_pointer(self): 

762 return self._first_pointer_parent[self._first_pointer_index] 

763 

764 @first_pointer.setter 

765 def first_pointer(self, value): 

766 self._first_pointer_parent[self._first_pointer_index] = value 

767 

768 @property 

769 def second_pointer(self): 

770 return self._second_pointer_parent[self._second_pointer_index] 

771 

772 @second_pointer.setter 

773 def second_pointer(self, value): 

774 self._second_pointer_parent[self._second_pointer_index] = value 

775 

776 def __repr__(self): 

777 built = [self._build(i) for i in self._lines] 

778 return "_LeftRightArgs(lines=%s, higher=%s)" % ( 

779 built, 

780 self.higher, 

781 ) 

782 

783 def transpose(self): 

784 self._first_pointer_parent, self._second_pointer_parent = self._second_pointer_parent, self._first_pointer_parent 

785 self._first_pointer_index, self._second_pointer_index = self._second_pointer_index, self._first_pointer_index 

786 self._first_line_index, self._second_line_index = self._second_line_index, self._first_line_index 

787 return self 

788 

789 @staticmethod 

790 def _build(expr): 

791 if isinstance(expr, ExprBuilder): 

792 return expr.build() 

793 if isinstance(expr, list): 

794 if len(expr) == 1: 

795 return expr[0] 

796 else: 

797 return expr[0](*[_LeftRightArgs._build(i) for i in expr[1]]) 

798 else: 

799 return expr 

800 

801 def build(self): 

802 data = [self._build(i) for i in self._lines] 

803 if self.higher != 1: 

804 data += [self._build(self.higher)] 

805 data = list(data)