__init__.py 8.1 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269270271
  1. r"""
  2. N-dim array module for SymPy.
  3. Four classes are provided to handle N-dim arrays, given by the combinations
  4. dense/sparse (i.e. whether to store all elements or only the non-zero ones in
  5. memory) and mutable/immutable (immutable classes are SymPy objects, but cannot
  6. change after they have been created).
  7. Examples
  8. ========
  9. The following examples show the usage of ``Array``. This is an abbreviation for
  10. ``ImmutableDenseNDimArray``, that is an immutable and dense N-dim array, the
  11. other classes are analogous. For mutable classes it is also possible to change
  12. element values after the object has been constructed.
  13. Array construction can detect the shape of nested lists and tuples:
  14. >>> from sympy import Array
  15. >>> a1 = Array([[1, 2], [3, 4], [5, 6]])
  16. >>> a1
  17. [[1, 2], [3, 4], [5, 6]]
  18. >>> a1.shape
  19. (3, 2)
  20. >>> a1.rank()
  21. 2
  22. >>> from sympy.abc import x, y, z
  23. >>> a2 = Array([[[x, y], [z, x*z]], [[1, x*y], [1/x, x/y]]])
  24. >>> a2
  25. [[[x, y], [z, x*z]], [[1, x*y], [1/x, x/y]]]
  26. >>> a2.shape
  27. (2, 2, 2)
  28. >>> a2.rank()
  29. 3
  30. Otherwise one could pass a 1-dim array followed by a shape tuple:
  31. >>> m1 = Array(range(12), (3, 4))
  32. >>> m1
  33. [[0, 1, 2, 3], [4, 5, 6, 7], [8, 9, 10, 11]]
  34. >>> m2 = Array(range(12), (3, 2, 2))
  35. >>> m2
  36. [[[0, 1], [2, 3]], [[4, 5], [6, 7]], [[8, 9], [10, 11]]]
  37. >>> m2[1,1,1]
  38. 7
  39. >>> m2.reshape(4, 3)
  40. [[0, 1, 2], [3, 4, 5], [6, 7, 8], [9, 10, 11]]
  41. Slice support:
  42. >>> m2[:, 1, 1]
  43. [3, 7, 11]
  44. Elementwise derivative:
  45. >>> from sympy.abc import x, y, z
  46. >>> m3 = Array([x**3, x*y, z])
  47. >>> m3.diff(x)
  48. [3*x**2, y, 0]
  49. >>> m3.diff(z)
  50. [0, 0, 1]
  51. Multiplication with other SymPy expressions is applied elementwisely:
  52. >>> (1+x)*m3
  53. [x**3*(x + 1), x*y*(x + 1), z*(x + 1)]
  54. To apply a function to each element of the N-dim array, use ``applyfunc``:
  55. >>> m3.applyfunc(lambda x: x/2)
  56. [x**3/2, x*y/2, z/2]
  57. N-dim arrays can be converted to nested lists by the ``tolist()`` method:
  58. >>> m2.tolist()
  59. [[[0, 1], [2, 3]], [[4, 5], [6, 7]], [[8, 9], [10, 11]]]
  60. >>> isinstance(m2.tolist(), list)
  61. True
  62. If the rank is 2, it is possible to convert them to matrices with ``tomatrix()``:
  63. >>> m1.tomatrix()
  64. Matrix([
  65. [0, 1, 2, 3],
  66. [4, 5, 6, 7],
  67. [8, 9, 10, 11]])
  68. Products and contractions
  69. -------------------------
  70. Tensor product between arrays `A_{i_1,\ldots,i_n}` and `B_{j_1,\ldots,j_m}`
  71. creates the combined array `P = A \otimes B` defined as
  72. `P_{i_1,\ldots,i_n,j_1,\ldots,j_m} := A_{i_1,\ldots,i_n}\cdot B_{j_1,\ldots,j_m}.`
  73. It is available through ``tensorproduct(...)``:
  74. >>> from sympy import Array, tensorproduct
  75. >>> from sympy.abc import x,y,z,t
  76. >>> A = Array([x, y, z, t])
  77. >>> B = Array([1, 2, 3, 4])
  78. >>> tensorproduct(A, B)
  79. [[x, 2*x, 3*x, 4*x], [y, 2*y, 3*y, 4*y], [z, 2*z, 3*z, 4*z], [t, 2*t, 3*t, 4*t]]
  80. In case you don't want to evaluate the tensor product immediately, you can use
  81. ``ArrayTensorProduct``, which creates an unevaluated tensor product expression:
  82. >>> from sympy.tensor.array.expressions import ArrayTensorProduct
  83. >>> ArrayTensorProduct(A, B)
  84. ArrayTensorProduct([x, y, z, t], [1, 2, 3, 4])
  85. Calling ``.as_explicit()`` on ``ArrayTensorProduct`` is equivalent to just calling
  86. ``tensorproduct(...)``:
  87. >>> ArrayTensorProduct(A, B).as_explicit()
  88. [[x, 2*x, 3*x, 4*x], [y, 2*y, 3*y, 4*y], [z, 2*z, 3*z, 4*z], [t, 2*t, 3*t, 4*t]]
  89. Tensor product between a rank-1 array and a matrix creates a rank-3 array:
  90. >>> from sympy import eye
  91. >>> p1 = tensorproduct(A, eye(4))
  92. >>> p1
  93. [[[x, 0, 0, 0], [0, x, 0, 0], [0, 0, x, 0], [0, 0, 0, x]], [[y, 0, 0, 0], [0, y, 0, 0], [0, 0, y, 0], [0, 0, 0, y]], [[z, 0, 0, 0], [0, z, 0, 0], [0, 0, z, 0], [0, 0, 0, z]], [[t, 0, 0, 0], [0, t, 0, 0], [0, 0, t, 0], [0, 0, 0, t]]]
  94. Now, to get back `A_0 \otimes \mathbf{1}` one can access `p_{0,m,n}` by slicing:
  95. >>> p1[0,:,:]
  96. [[x, 0, 0, 0], [0, x, 0, 0], [0, 0, x, 0], [0, 0, 0, x]]
  97. Tensor contraction sums over the specified axes, for example contracting
  98. positions `a` and `b` means
  99. `A_{i_1,\ldots,i_a,\ldots,i_b,\ldots,i_n} \implies \sum_k A_{i_1,\ldots,k,\ldots,k,\ldots,i_n}`
  100. Remember that Python indexing is zero starting, to contract the a-th and b-th
  101. axes it is therefore necessary to specify `a-1` and `b-1`
  102. >>> from sympy import tensorcontraction
  103. >>> C = Array([[x, y], [z, t]])
  104. The matrix trace is equivalent to the contraction of a rank-2 array:
  105. `A_{m,n} \implies \sum_k A_{k,k}`
  106. >>> tensorcontraction(C, (0, 1))
  107. t + x
  108. To create an expression representing a tensor contraction that does not get
  109. evaluated immediately, use ``ArrayContraction``, which is equivalent to
  110. ``tensorcontraction(...)`` if it is followed by ``.as_explicit()``:
  111. >>> from sympy.tensor.array.expressions import ArrayContraction
  112. >>> ArrayContraction(C, (0, 1))
  113. ArrayContraction([[x, y], [z, t]], (0, 1))
  114. >>> ArrayContraction(C, (0, 1)).as_explicit()
  115. t + x
  116. Matrix product is equivalent to a tensor product of two rank-2 arrays, followed
  117. by a contraction of the 2nd and 3rd axes (in Python indexing axes number 1, 2).
  118. `A_{m,n}\cdot B_{i,j} \implies \sum_k A_{m, k}\cdot B_{k, j}`
  119. >>> D = Array([[2, 1], [0, -1]])
  120. >>> tensorcontraction(tensorproduct(C, D), (1, 2))
  121. [[2*x, x - y], [2*z, -t + z]]
  122. One may verify that the matrix product is equivalent:
  123. >>> from sympy import Matrix
  124. >>> Matrix([[x, y], [z, t]])*Matrix([[2, 1], [0, -1]])
  125. Matrix([
  126. [2*x, x - y],
  127. [2*z, -t + z]])
  128. or equivalently
  129. >>> C.tomatrix()*D.tomatrix()
  130. Matrix([
  131. [2*x, x - y],
  132. [2*z, -t + z]])
  133. Diagonal operator
  134. -----------------
  135. The ``tensordiagonal`` function acts in a similar manner as ``tensorcontraction``,
  136. but the joined indices are not summed over, for example diagonalizing
  137. positions `a` and `b` means
  138. `A_{i_1,\ldots,i_a,\ldots,i_b,\ldots,i_n} \implies A_{i_1,\ldots,k,\ldots,k,\ldots,i_n}
  139. \implies \tilde{A}_{i_1,\ldots,i_{a-1},i_{a+1},\ldots,i_{b-1},i_{b+1},\ldots,i_n,k}`
  140. where `\tilde{A}` is the array equivalent to the diagonal of `A` at positions
  141. `a` and `b` moved to the last index slot.
  142. Compare the difference between contraction and diagonal operators:
  143. >>> from sympy import tensordiagonal
  144. >>> from sympy.abc import a, b, c, d
  145. >>> m = Matrix([[a, b], [c, d]])
  146. >>> tensorcontraction(m, [0, 1])
  147. a + d
  148. >>> tensordiagonal(m, [0, 1])
  149. [a, d]
  150. In short, no summation occurs with ``tensordiagonal``.
  151. Derivatives by array
  152. --------------------
  153. The usual derivative operation may be extended to support derivation with
  154. respect to arrays, provided that all elements in the that array are symbols or
  155. expressions suitable for derivations.
  156. The definition of a derivative by an array is as follows: given the array
  157. `A_{i_1, \ldots, i_N}` and the array `X_{j_1, \ldots, j_M}`
  158. the derivative of arrays will return a new array `B` defined by
  159. `B_{j_1,\ldots,j_M,i_1,\ldots,i_N} := \frac{\partial A_{i_1,\ldots,i_N}}{\partial X_{j_1,\ldots,j_M}}`
  160. The function ``derive_by_array`` performs such an operation:
  161. >>> from sympy import derive_by_array
  162. >>> from sympy.abc import x, y, z, t
  163. >>> from sympy import sin, exp
  164. With scalars, it behaves exactly as the ordinary derivative:
  165. >>> derive_by_array(sin(x*y), x)
  166. y*cos(x*y)
  167. Scalar derived by an array basis:
  168. >>> derive_by_array(sin(x*y), [x, y, z])
  169. [y*cos(x*y), x*cos(x*y), 0]
  170. Deriving array by an array basis: `B^{nm} := \frac{\partial A^m}{\partial x^n}`
  171. >>> basis = [x, y, z]
  172. >>> ax = derive_by_array([exp(x), sin(y*z), t], basis)
  173. >>> ax
  174. [[exp(x), 0, 0], [0, z*cos(y*z), 0], [0, y*cos(y*z), 0]]
  175. Contraction of the resulting array: `\sum_m \frac{\partial A^m}{\partial x^m}`
  176. >>> tensorcontraction(ax, (0, 1))
  177. z*cos(y*z) + exp(x)
  178. """
  179. from .dense_ndim_array import MutableDenseNDimArray, ImmutableDenseNDimArray, DenseNDimArray
  180. from .sparse_ndim_array import MutableSparseNDimArray, ImmutableSparseNDimArray, SparseNDimArray
  181. from .ndim_array import NDimArray, ArrayKind
  182. from .arrayop import tensorproduct, tensorcontraction, tensordiagonal, derive_by_array, permutedims
  183. from .array_comprehension import ArrayComprehension, ArrayComprehensionMap
  184. Array = ImmutableDenseNDimArray
  185. __all__ = [
  186. 'MutableDenseNDimArray', 'ImmutableDenseNDimArray', 'DenseNDimArray',
  187. 'MutableSparseNDimArray', 'ImmutableSparseNDimArray', 'SparseNDimArray',
  188. 'NDimArray', 'ArrayKind',
  189. 'tensorproduct', 'tensorcontraction', 'tensordiagonal', 'derive_by_array',
  190. 'permutedims', 'ArrayComprehension', 'ArrayComprehensionMap',
  191. 'Array',
  192. ]