import random
from sympy.combinatorics import Permutation
from sympy.combinatorics.permutations import _af_invert
from sympy.testing.pytest import raises
from sympy.core.function import diff
from sympy.core.symbol import symbols
from sympy.functions.elementary.complexes import (adjoint, conjugate, transpose)
from sympy.functions.elementary.exponential import (exp, log)
from sympy.functions.elementary.trigonometric import (cos, sin)
from sympy.tensor.array import Array, ImmutableDenseNDimArray, ImmutableSparseNDimArray, MutableSparseNDimArray
from sympy.tensor.array.arrayop import tensorproduct, tensorcontraction, derive_by_array, permutedims, Flatten, \
tensordiagonal
def test_import_NDimArray():
from sympy.tensor.array import NDimArray
del NDimArray
def test_tensorproduct():
x,y,z,t = symbols('x y z t')
from sympy.abc import a,b,c,d
assert tensorproduct() == 1
assert tensorproduct([x]) == Array([x])
assert tensorproduct([x], [y]) == Array([[x*y]])
assert tensorproduct([x], [y], [z]) == Array([[[x*y*z]]])
assert tensorproduct([x], [y], [z], [t]) == Array([[[[x*y*z*t]]]])
assert tensorproduct(x) == x
assert tensorproduct(x, y) == x*y
assert tensorproduct(x, y, z) == x*y*z
assert tensorproduct(x, y, z, t) == x*y*z*t
for ArrayType in [ImmutableDenseNDimArray, ImmutableSparseNDimArray]:
A = ArrayType([x, y])
B = ArrayType([1, 2, 3])
C = ArrayType([a, b, c, d])
assert tensorproduct(A, B, C) == ArrayType([[[a*x, b*x, c*x, d*x], [2*a*x, 2*b*x, 2*c*x, 2*d*x], [3*a*x, 3*b*x, 3*c*x, 3*d*x]],
[[a*y, b*y, c*y, d*y], [2*a*y, 2*b*y, 2*c*y, 2*d*y], [3*a*y, 3*b*y, 3*c*y, 3*d*y]]])
assert tensorproduct([x, y], [1, 2, 3]) == tensorproduct(A, B)
assert tensorproduct(A, 2) == ArrayType([2*x, 2*y])
assert tensorproduct(A, [2]) == ArrayType([[2*x], [2*y]])
assert tensorproduct([2], A) == ArrayType([[2*x, 2*y]])
assert tensorproduct(a, A) == ArrayType([a*x, a*y])
assert tensorproduct(a, A, B) == ArrayType([[a*x, 2*a*x, 3*a*x], [a*y, 2*a*y, 3*a*y]])
assert tensorproduct(A, B, a) == ArrayType([[a*x, 2*a*x, 3*a*x], [a*y, 2*a*y, 3*a*y]])
assert tensorproduct(B, a, A) == ArrayType([[a*x, a*y], [2*a*x, 2*a*y], [3*a*x, 3*a*y]])
# tests for large scale sparse array
for SparseArrayType in [ImmutableSparseNDimArray, MutableSparseNDimArray]:
a = SparseArrayType({1:2, 3:4},(1000, 2000))
b = SparseArrayType({1:2, 3:4},(1000, 2000))
assert tensorproduct(a, b) == ImmutableSparseNDimArray({2000001: 4, 2000003: 8, 6000001: 8, 6000003: 16}, (1000, 2000, 1000, 2000))
def test_tensorcontraction():
from sympy.abc import a,b,c,d,e,f,g,h,i,j,k,l,m,n,o,p,q,r,s,t,u,v,w,x
B = Array(range(18), (2, 3, 3))
assert tensorcontraction(B, (1, 2)) == Array([12, 39])
C1 = Array([a, b, c, d, e, f, g, h, i, j, k, l, m, n, o, p, q, r, s, t, u, v, w, x], (2, 3, 2, 2))
assert tensorcontraction(C1, (0, 2)) == Array([[a + o, b + p], [e + s, f + t], [i + w, j + x]])
assert tensorcontraction(C1, (0, 2, 3)) == Array([a + p, e + t, i + x])
assert tensorcontraction(C1, (2, 3)) == Array([[a + d, e + h, i + l], [m + p, q + t, u + x]])
def test_derivative_by_array():
from sympy.abc import i, j, t, x, y, z
bexpr = x*y**2*exp(z)*log(t)
sexpr = sin(bexpr)
cexpr = cos(bexpr)
a = Array([sexpr])
assert derive_by_array(sexpr, t) == x*y**2*exp(z)*cos(x*y**2*exp(z)*log(t))/t
assert derive_by_array(sexpr, [x, y, z]) == Array([bexpr/x*cexpr, 2*y*bexpr/y**2*cexpr, bexpr*cexpr])
assert derive_by_array(a, [x, y, z]) == Array([[bexpr/x*cexpr], [2*y*bexpr/y**2*cexpr], [bexpr*cexpr]])
assert derive_by_array(sexpr, [[x, y], [z, t]]) == Array([[bexpr/x*cexpr, 2*y*bexpr/y**2*cexpr], [bexpr*cexpr, bexpr/log(t)/t*cexpr]])
assert derive_by_array(a, [[x, y], [z, t]]) == Array([[[bexpr/x*cexpr], [2*y*bexpr/y**2*cexpr]], [[bexpr*cexpr], [bexpr/log(t)/t*cexpr]]])
assert derive_by_array([[x, y], [z, t]], [x, y]) == Array([[[1, 0], [0, 0]], [[0, 1], [0, 0]]])
assert derive_by_array([[x, y], [z, t]], [[x, y], [z, t]]) == Array([[[[1, 0], [0, 0]], [[0, 1], [0, 0]]],
[[[0, 0], [1, 0]], [[0, 0], [0, 1]]]])
assert diff(sexpr, t) == x*y**2*exp(z)*cos(x*y**2*exp(z)*log(t))/t
assert diff(sexpr, Array([x, y, z])) == Array([bexpr/x*cexpr, 2*y*bexpr/y**2*cexpr, bexpr*cexpr])
assert diff(a, Array([x, y, z])) == Array([[bexpr/x*cexpr], [2*y*bexpr/y**2*cexpr], [bexpr*cexpr]])
assert diff(sexpr, Array([[x, y], [z, t]])) == Array([[bexpr/x*cexpr, 2*y*bexpr/y**2*cexpr], [bexpr*cexpr, bexpr/log(t)/t*cexpr]])
assert diff(a, Array([[x, y], [z, t]])) == Array([[[bexpr/x*cexpr], [2*y*bexpr/y**2*cexpr]], [[bexpr*cexpr], [bexpr/log(t)/t*cexpr]]])
assert diff(Array([[x, y], [z, t]]), Array([x, y])) == Array([[[1, 0], [0, 0]], [[0, 1], [0, 0]]])
assert diff(Array([[x, y], [z, t]]), Array([[x, y], [z, t]])) == Array([[[[1, 0], [0, 0]], [[0, 1], [0, 0]]],
[[[0, 0], [1, 0]], [[0, 0], [0, 1]]]])
# test for large scale sparse array
for SparseArrayType in [ImmutableSparseNDimArray, MutableSparseNDimArray]:
b = MutableSparseNDimArray({0:i, 1:j}, (10000, 20000))
assert derive_by_array(b, i) == ImmutableSparseNDimArray({0: 1}, (10000, 20000))
assert derive_by_array(b, (i, j)) == ImmutableSparseNDimArray({0: 1, 200000001: 1}, (2, 10000, 20000))
#https://github.com/sympy/sympy/issues/20655
U = Array([x, y, z])
E = 2
assert derive_by_array(E, U) == ImmutableDenseNDimArray([0, 0, 0])
def test_issue_emerged_while_discussing_10972():
ua = Array([-1,0])
Fa = Array([[0, 1], [-1, 0]])
po = tensorproduct(Fa, ua, Fa, ua)
assert tensorcontraction(po, (1, 2), (4, 5)) == Array([[0, 0], [0, 1]])
sa = symbols('a0:144')
po = Array(sa, [2, 2, 3, 3, 2, 2])
assert tensorcontraction(po, (0, 1), (2, 3), (4, 5)) == sa[0] + sa[108] + sa[111] + sa[124] + sa[127] + sa[140] + sa[143] + sa[16] + sa[19] + sa[3] + sa[32] + sa[35]
assert tensorcontraction(po, (0, 1, 4, 5), (2, 3)) == sa[0] + sa[111] + sa[127] + sa[143] + sa[16] + sa[32]
assert tensorcontraction(po, (0, 1), (4, 5)) == Array([[sa[0] + sa[108] + sa[111] + sa[3], sa[112] + sa[115] + sa[4] + sa[7],
sa[11] + sa[116] + sa[119] + sa[8]], [sa[12] + sa[120] + sa[123] + sa[15],
sa[124] + sa[127] + sa[16] + sa[19], sa[128] + sa[131] + sa[20] + sa[23]],
[sa[132] + sa[135] + sa[24] + sa[27], sa[136] + sa[139] + sa[28] + sa[31],
sa[140] + sa[143] + sa[32] + sa[35]]])
assert tensorcontraction(po, (0, 1), (2, 3)) == Array([[sa[0] + sa[108] + sa[124] + sa[140] + sa[16] + sa[32], sa[1] + sa[109] + sa[125] + sa[141] + sa[17] + sa[33]],
[sa[110] + sa[126] + sa[142] + sa[18] + sa[2] + sa[34], sa[111] + sa[127] + sa[143] + sa[19] + sa[3] + sa[35]]])
def test_array_permutedims():
sa = symbols('a0:144')
for ArrayType in [ImmutableDenseNDimArray, ImmutableSparseNDimArray]:
m1 = ArrayType(sa[:6], (2, 3))
assert permutedims(m1, (1, 0)) == transpose(m1)
assert m1.tomatrix().T == permutedims(m1, (1, 0)).tomatrix()
assert m1.tomatrix().T == transpose(m1).tomatrix()
assert m1.tomatrix().C == conjugate(m1).tomatrix()
assert m1.tomatrix().H == adjoint(m1).tomatrix()
assert m1.tomatrix().T == m1.transpose().tomatrix()
assert m1.tomatrix().C == m1.conjugate().tomatrix()
assert m1.tomatrix().H == m1.adjoint().tomatrix()
raises(ValueError, lambda: permutedims(m1, (0,)))
raises(ValueError, lambda: permutedims(m1, (0, 0)))
raises(ValueError, lambda: permutedims(m1, (1, 2, 0)))
# Some tests with random arrays:
dims = 6
shape = [random.randint(1,5) for i in range(dims)]
elems = [random.random() for i in range(tensorproduct(*shape))]
ra = ArrayType(elems, shape)
perm = list(range(dims))
# Randomize the permutation:
random.shuffle(perm)
# Test inverse permutation:
assert permutedims(permutedims(ra, perm), _af_invert(perm)) == ra
# Test that permuted shape corresponds to action by `Permutation`:
assert permutedims(ra, perm).shape == tuple(Permutation(perm)(shape))
z = ArrayType.zeros(4,5,6,7)
assert permutedims(z, (2, 3, 1, 0)).shape == (6, 7, 5, 4)
assert permutedims(z, [2, 3, 1, 0]).shape == (6, 7, 5, 4)
assert permutedims(z, Permutation([2, 3, 1, 0])).shape == (6, 7, 5, 4)
po = ArrayType(sa, [2, 2, 3, 3, 2, 2])
raises(ValueError, lambda: permutedims(po, (1, 1)))
raises(ValueError, lambda: po.transpose())
raises(ValueError, lambda: po.adjoint())
assert permutedims(po, reversed(range(po.rank()))) == ArrayType(
[[[[[[sa[0], sa[72]], [sa[36], sa[108]]], [[sa[12], sa[84]], [sa[48], sa[120]]], [[sa[24],
sa[96]], [sa[60], sa[132]]]],
[[[sa[4], sa[76]], [sa[40], sa[112]]], [[sa[16],
sa[88]], [sa[52], sa[124]]],
[[sa[28], sa[100]], [sa[64], sa[136]]]],
[[[sa[8],
sa[80]], [sa[44], sa[116]]], [[sa[20], sa[92]], [sa[56], sa[128]]], [[sa[32],
sa[104]], [sa[68], sa[140]]]]],
[[[[sa[2], sa[74]], [sa[38], sa[110]]], [[sa[14],
sa[86]], [sa[50], sa[122]]], [[sa[26], sa[98]], [sa[62], sa[134]]]],
[[[sa[6],
sa[78]], [sa[42], sa[114]]], [[sa[18], sa[90]], [sa[54], sa[126]]], [[sa[30],
sa[102]], [sa[66], sa[138]]]],
[[[sa[10], sa[82]], [sa[46], sa[118]]], [[sa[22],
sa[94]], [sa[58], sa[130]]],
[[sa[34], sa[106]], [sa[70], sa[142]]]]]],
[[[[[sa[1],
sa[73]], [sa[37], sa[109]]], [[sa[13], sa[85]], [sa[49], sa[121]]], [[sa[25],
sa[97]], [sa[61], sa[133]]]],
[[[sa[5], sa[77]], [sa[41], sa[113]]], [[sa[17],
sa[89]], [sa[53], sa[125]]],
[[sa[29], sa[101]], [sa[65], sa[137]]]],
[[[sa[9],
sa[81]], [sa[45], sa[117]]], [[sa[21], sa[93]], [sa[57], sa[129]]], [[sa[33],
sa[105]], [sa[69], sa[141]]]]],
[[[[sa[3], sa[75]], [sa[39], sa[111]]], [[sa[15],
sa[87]], [sa[51], sa[123]]], [[sa[27], sa[99]], [sa[63], sa[135]]]],
[[[sa[7],
sa[79]], [sa[43], sa[115]]], [[sa[19], sa[91]], [sa[55], sa[127]]], [[sa[31],
sa[103]], [sa[67], sa[139]]]],
[[[sa[11], sa[83]], [sa[47], sa[119]]], [[sa[23],
sa[95]], [sa[59], sa[131]]],
[[sa[35], sa[107]], [sa[71], sa[143]]]]]]])
assert permutedims(po, (1, 0, 2, 3, 4, 5)) == ArrayType(
[[[[[[sa[0], sa[1]], [sa[2], sa[3]]], [[sa[4], sa[5]], [sa[6], sa[7]]], [[sa[8], sa[9]], [sa[10],
sa[11]]]],
[[[sa[12], sa[13]], [sa[14], sa[15]]], [[sa[16], sa[17]], [sa[18],
sa[19]]], [[sa[20], sa[21]], [sa[22], sa[23]]]],
[[[sa[24], sa[25]], [sa[26],
sa[27]]], [[sa[28], sa[29]], [sa[30], sa[31]]], [[sa[32], sa[33]], [sa[34],
sa[35]]]]],
[[[[sa[72], sa[73]], [sa[74], sa[75]]], [[sa[76], sa[77]], [sa[78],
sa[79]]], [[sa[80], sa[81]], [sa[82], sa[83]]]],
[[[sa[84], sa[85]], [sa[86],
sa[87]]], [[sa[88], sa[89]], [sa[90], sa[91]]], [[sa[92], sa[93]], [sa[94],
sa[95]]]],
[[[sa[96], sa[97]], [sa[98], sa[99]]], [[sa[100], sa[101]], [sa[102],
sa[103]]],
[[sa[104], sa[105]], [sa[106], sa[107]]]]]], [[[[[sa[36], sa[37]], [sa[38],
sa[39]]],
[[sa[40], sa[41]], [sa[42], sa[43]]],
[[sa[44], sa[45]], [sa[46],
sa[47]]]],
[[[sa[48], sa[49]], [sa[50], sa[51]]],
[[sa[52], sa[53]], [sa[54],
sa[55]]],
[[sa[56], sa[57]], [sa[58], sa[59]]]],
[[[sa[60], sa[61]], [sa[62],
sa[63]]],
[[sa[64], sa[65]], [sa[66], sa[67]]],
[[sa[68], sa[69]], [sa[70],
sa[71]]]]], [
[[[sa[108], sa[109]], [sa[110], sa[111]]],
[[sa[112], sa[113]], [sa[114],
sa[115]]],
[[sa[116], sa[117]], [sa[118], sa[119]]]],
[[[sa[120], sa[121]], [sa[122],
sa[123]]],
[[sa[124], sa[125]], [sa[126], sa[127]]],
[[sa[128], sa[129]], [sa[130],
sa[131]]]],
[[[sa[132], sa[133]], [sa[134], sa[135]]],
[[sa[136], sa[137]], [sa[138],
sa[139]]],
[[sa[140], sa[141]], [sa[142], sa[143]]]]]]])
assert permutedims(po, (0, 2, 1, 4, 3, 5)) == ArrayType(
[[[[[[sa[0], sa[1]], [sa[4], sa[5]], [sa[8], sa[9]]], [[sa[2], sa[3]], [sa[6], sa[7]], [sa[10],
sa[11]]]],
[[[sa[36], sa[37]], [sa[40], sa[41]], [sa[44], sa[45]]], [[sa[38],
sa[39]], [sa[42], sa[43]], [sa[46], sa[47]]]]],
[[[[sa[12], sa[13]], [sa[16],
sa[17]], [sa[20], sa[21]]], [[sa[14], sa[15]], [sa[18], sa[19]], [sa[22],
sa[23]]]],
[[[sa[48], sa[49]], [sa[52], sa[53]], [sa[56], sa[57]]], [[sa[50],
sa[51]], [sa[54], sa[55]], [sa[58], sa[59]]]]],
[[[[sa[24], sa[25]], [sa[28],
sa[29]], [sa[32], sa[33]]], [[sa[26], sa[27]], [sa[30], sa[31]], [sa[34],
sa[35]]]],
[[[sa[60], sa[61]], [sa[64], sa[65]], [sa[68], sa[69]]], [[sa[62],
sa[63]], [sa[66], sa[67]], [sa[70], sa[71]]]]]],
[[[[[sa[72], sa[73]], [sa[76],
sa[77]], [sa[80], sa[81]]], [[sa[74], sa[75]], [sa[78], sa[79]], [sa[82],
sa[83]]]],
[[[sa[108], sa[109]], [sa[112], sa[113]], [sa[116], sa[117]]], [[sa[110],
sa[111]], [sa[114], sa[115]],
[sa[118], sa[119]]]]],
[[[[sa[84], sa[85]], [sa[88],
sa[89]], [sa[92], sa[93]]], [[sa[86], sa[87]], [sa[90], sa[91]], [sa[94],
sa[95]]]],
[[[sa[120], sa[121]], [sa[124], sa[125]], [sa[128], sa[129]]], [[sa[122],
sa[123]], [sa[126], sa[127]],
[sa[130], sa[131]]]]],
[[[[sa[96], sa[97]], [sa[100],
sa[101]], [sa[104], sa[105]]], [[sa[98], sa[99]], [sa[102], sa[103]], [sa[106],
sa[107]]]],
[[[sa[132], sa[133]], [sa[136], sa[137]], [sa[140], sa[141]]], [[sa[134],
sa[135]], [sa[138], sa[139]],
[sa[142], sa[143]]]]]]])
po2 = po.reshape(4, 9, 2, 2)
assert po2 == ArrayType([[[[sa[0], sa[1]], [sa[2], sa[3]]], [[sa[4], sa[5]], [sa[6], sa[7]]], [[sa[8], sa[9]], [sa[10], sa[11]]], [[sa[12], sa[13]], [sa[14], sa[15]]], [[sa[16], sa[17]], [sa[18], sa[19]]], [[sa[20], sa[21]], [sa[22], sa[23]]], [[sa[24], sa[25]], [sa[26], sa[27]]], [[sa[28], sa[29]], [sa[30], sa[31]]], [[sa[32], sa[33]], [sa[34], sa[35]]]], [[[sa[36], sa[37]], [sa[38], sa[39]]], [[sa[40], sa[41]], [sa[42], sa[43]]], [[sa[44], sa[45]], [sa[46], sa[47]]], [[sa[48], sa[49]], [sa[50], sa[51]]], [[sa[52], sa[53]], [sa[54], sa[55]]], [[sa[56], sa[57]], [sa[58], sa[59]]], [[sa[60], sa[61]], [sa[62], sa[63]]], [[sa[64], sa[65]], [sa[66], sa[67]]], [[sa[68], sa[69]], [sa[70], sa[71]]]], [[[sa[72], sa[73]], [sa[74], sa[75]]], [[sa[76], sa[77]], [sa[78], sa[79]]], [[sa[80], sa[81]], [sa[82], sa[83]]], [[sa[84], sa[85]], [sa[86], sa[87]]], [[sa[88], sa[89]], [sa[90], sa[91]]], [[sa[92], sa[93]], [sa[94], sa[95]]], [[sa[96], sa[97]], [sa[98], sa[99]]], [[sa[100], sa[101]], [sa[102], sa[103]]], [[sa[104], sa[105]], [sa[106], sa[107]]]], [[[sa[108], sa[109]], [sa[110], sa[111]]], [[sa[112], sa[113]], [sa[114], sa[115]]], [[sa[116], sa[117]], [sa[118], sa[119]]], [[sa[120], sa[121]], [sa[122], sa[123]]], [[sa[124], sa[125]], [sa[126], sa[127]]], [[sa[128], sa[129]], [sa[130], sa[131]]], [[sa[132], sa[133]], [sa[134], sa[135]]], [[sa[136], sa[137]], [sa[138], sa[139]]], [[sa[140], sa[141]], [sa[142], sa[143]]]]])
assert permutedims(po2, (3, 2, 0, 1)) == ArrayType([[[[sa[0], sa[4], sa[8], sa[12], sa[16], sa[20], sa[24], sa[28], sa[32]], [sa[36], sa[40], sa[44], sa[48], sa[52], sa[56], sa[60], sa[64], sa[68]], [sa[72], sa[76], sa[80], sa[84], sa[88], sa[92], sa[96], sa[100], sa[104]], [sa[108], sa[112], sa[116], sa[120], sa[124], sa[128], sa[132], sa[136], sa[140]]], [[sa[2], sa[6], sa[10], sa[14], sa[18], sa[22], sa[26], sa[30], sa[34]], [sa[38], sa[42], sa[46], sa[50], sa[54], sa[58], sa[62], sa[66], sa[70]], [sa[74], sa[78], sa[82], sa[86], sa[90], sa[94], sa[98], sa[102], sa[106]], [sa[110], sa[114], sa[118], sa[122], sa[126], sa[130], sa[134], sa[138], sa[142]]]], [[[sa[1], sa[5], sa[9], sa[13], sa[17], sa[21], sa[25], sa[29], sa[33]], [sa[37], sa[41], sa[45], sa[49], sa[53], sa[57], sa[61], sa[65], sa[69]], [sa[73], sa[77], sa[81], sa[85], sa[89], sa[93], sa[97], sa[101], sa[105]], [sa[109], sa[113], sa[117], sa[121], sa[125], sa[129], sa[133], sa[137], sa[141]]], [[sa[3], sa[7], sa[11], sa[15], sa[19], sa[23], sa[27], sa[31], sa[35]], [sa[39], sa[43], sa[47], sa[51], sa[55], sa[59], sa[63], sa[67], sa[71]], [sa[75], sa[79], sa[83], sa[87], sa[91], sa[95], sa[99], sa[103], sa[107]], [sa[111], sa[115], sa[119], sa[123], sa[127], sa[131], sa[135], sa[139], sa[143]]]]])
# test for large scale sparse array
for SparseArrayType in [ImmutableSparseNDimArray, MutableSparseNDimArray]:
A = SparseArrayType({1:1, 10000:2}, (10000, 20000, 10000))
assert permutedims(A, (0, 1, 2)) == A
assert permutedims(A, (1, 0, 2)) == SparseArrayType({1: 1, 100000000: 2}, (20000, 10000, 10000))
B = SparseArrayType({1:1, 20000:2}, (10000, 20000))
assert B.transpose() == SparseArrayType({10000: 1, 1: 2}, (20000, 10000))
def test_flatten():
from sympy.matrices.dense import Matrix
for ArrayType in [ImmutableDenseNDimArray, ImmutableSparseNDimArray, Matrix]:
A = ArrayType(range(24)).reshape(4, 6)
assert [i for i in Flatten(A)] == [0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23]
for i, v in enumerate(Flatten(A)):
assert i == v
def test_tensordiagonal():
from sympy.matrices.dense import eye
expr = Array(range(9)).reshape(3, 3)
raises(ValueError, lambda: tensordiagonal(expr, [0], [1]))
raises(ValueError, lambda: tensordiagonal(expr, [0, 0]))
assert tensordiagonal(eye(3), [0, 1]) == Array([1, 1, 1])
assert tensordiagonal(expr, [0, 1]) == Array([0, 4, 8])
x, y, z = symbols("x y z")
expr2 = tensorproduct([x, y, z], expr)
assert tensordiagonal(expr2, [1, 2]) == Array([[0, 4*x, 8*x], [0, 4*y, 8*y], [0, 4*z, 8*z]])
assert tensordiagonal(expr2, [0, 1]) == Array([[0, 3*y, 6*z], [x, 4*y, 7*z], [2*x, 5*y, 8*z]])
assert tensordiagonal(expr2, [0, 1, 2]) == Array([0, 4*y, 8*z])
# assert tensordiagonal(expr2, [0]) == permutedims(expr2, [1, 2, 0])
# assert tensordiagonal(expr2, [1]) == permutedims(expr2, [0, 2, 1])
# assert tensordiagonal(expr2, [2]) == expr2
# assert tensordiagonal(expr2, [1], [2]) == expr2
# assert tensordiagonal(expr2, [0], [1]) == permutedims(expr2, [2, 0, 1])
a, b, c, X, Y, Z = symbols("a b c X Y Z")
expr3 = tensorproduct([x, y, z], [1, 2, 3], [a, b, c], [X, Y, Z])
assert tensordiagonal(expr3, [0, 1, 2, 3]) == Array([x*a*X, 2*y*b*Y, 3*z*c*Z])
assert tensordiagonal(expr3, [0, 1], [2, 3]) == tensorproduct([x, 2*y, 3*z], [a*X, b*Y, c*Z])
# assert tensordiagonal(expr3, [0], [1, 2], [3]) == tensorproduct([x, y, z], [a, 2*b, 3*c], [X, Y, Z])
assert tensordiagonal(tensordiagonal(expr3, [2, 3]), [0, 1]) == tensorproduct([a*X, b*Y, c*Z], [x, 2*y, 3*z])
raises(ValueError, lambda: tensordiagonal([[1, 2, 3], [4, 5, 6]], [0, 1]))
raises(ValueError, lambda: tensordiagonal(expr3.reshape(3, 3, 9), [1, 2]))