# python - style - 使用NumPy進行快速張量旋轉

## 前瞻性方法和解決方案代碼

``````np.einsum('ai,bj,ck,dl,abcd->ijkl', g, g, g, g, T)
``````

``````p1 = np.einsum('abcd, ai->bcdi', T, g)
``````

``````p2 = np.einsum('bcdi, bj->cdij', p1, g)
``````

``````p3 = np.einsum('cdij, ck->dijk', p2, g)

p4 = np.einsum('dijk, dl->ijkl', p3, g)
``````

``````p1 = np.tensordot(T,g,axes=((0),(0)))
p2 = np.tensordot(p1,g,axes=((0),(0)))
p3 = np.tensordot(p2,g,axes=((0),(0)))
out = np.tensordot(p3,g,axes=((0),(0)))
``````

## 運行時測試

``````def rotT_Philipp(T, g):  # @Philipp's soln
gg = np.outer(g, g)
gggg = np.outer(gg, gg).reshape(4 * g.shape)
axes = ((0, 2, 4, 6), (0, 1, 2, 3))
return np.tensordot(gggg, T, axes)

def rotT_Sven(T, g):    # @Sven Marnach's soln
Tprime = T
for i in range(4):
slices = [None] * 4
slices[i] = slice(None)
slices *= 2
Tprime = g[slices].T * Tprime
return Tprime.sum(-1).sum(-1).sum(-1).sum(-1)

def rotT_pv(T, g):     # @pv.'s soln
return np.einsum('ai,bj,ck,dl,abcd->ijkl', g, g, g, g, T)

def rotT_Divakar(T,g): # Posted in this post
p1 = np.tensordot(T,g,axes=((0),(0)))
p2 = np.tensordot(p1,g,axes=((0),(0)))
p3 = np.tensordot(p2,g,axes=((0),(0)))
p4 = np.tensordot(p3,g,axes=((0),(0)))
return p4
``````

``````In [304]: # Setup inputs
...: T = np.random.rand(3,3,3,3)
...: g = np.random.rand(3,3)
...:

In [305]: %timeit rotT(T, g)
...: %timeit rotT_pv(T, g)
...: %timeit rotT_Sven(T, g)
...: %timeit rotT_Philipp(T, g)
...: %timeit rotT_Divakar(T, g)
...:
100 loops, best of 3: 6.51 ms per loop
1000 loops, best of 3: 247 µs per loop
10000 loops, best of 3: 137 µs per loop
10000 loops, best of 3: 41.6 µs per loop
10000 loops, best of 3: 28.3 µs per loop

In [306]: 6510.0/28.3 # Speedup with the proposed soln over original code
Out[306]: 230.03533568904592
``````

``````In [307]: # Setup inputs
...: T = np.random.rand(5,5,5,5)
...: g = np.random.rand(5,5)
...:

In [308]: %timeit rotT(T, g)
...: %timeit rotT_pv(T, g)
...: %timeit rotT_Sven(T, g)
...: %timeit rotT_Philipp(T, g)
...: %timeit rotT_Divakar(T, g)
...:
100 loops, best of 3: 6.54 ms per loop
100 loops, best of 3: 7.17 ms per loop
100 loops, best of 3: 2.7 ms per loop
1000 loops, best of 3: 1.47 ms per loop
10000 loops, best of 3: 39.9 µs per loop
``````

``````#!/usr/bin/env python

import numpy as np

def rotT(T, g):
Tprime = np.zeros((3,3,3,3))
for i in range(3):
for j in range(3):
for k in range(3):
for l in range(3):
for ii in range(3):
for jj in range(3):
for kk in range(3):
for ll in range(3):
gg = g[ii,i]*g[jj,j]*g[kk,k]*g[ll,l]
Tprime[i,j,k,l] = Tprime[i,j,k,l] + \
gg*T[ii,jj,kk,ll]
return Tprime

if __name__ == "__main__":

T = np.array([[[[  4.66533067e+01,  5.84985000e-02, -5.37671310e-01],
[  5.84985000e-02,  1.56722231e+01,  2.32831900e-02],
[ -5.37671310e-01,  2.32831900e-02,  1.33399259e+01]],
[[  4.60051700e-02,  1.54658176e+01,  2.19568200e-02],
[  1.54658176e+01, -5.18223500e-02, -1.52814920e-01],
[  2.19568200e-02, -1.52814920e-01, -2.43874100e-02]],
[[ -5.35577630e-01,  1.95558600e-02,  1.31108757e+01],
[  1.95558600e-02, -1.51342210e-01, -6.67615000e-03],
[  1.31108757e+01, -6.67615000e-03,  6.90486240e-01]]],
[[[  4.60051700e-02,  1.54658176e+01,  2.19568200e-02],
[  1.54658176e+01, -5.18223500e-02, -1.52814920e-01],
[  2.19568200e-02, -1.52814920e-01, -2.43874100e-02]],
[[  1.57414726e+01, -3.86167500e-02, -1.55971950e-01],
[ -3.86167500e-02,  4.65601977e+01, -3.57741000e-02],
[ -1.55971950e-01, -3.57741000e-02,  1.34215636e+01]],
[[  2.58256300e-02, -1.49072770e-01, -7.38843000e-03],
[ -1.49072770e-01, -3.63410500e-02,  1.32039847e+01],
[ -7.38843000e-03,  1.32039847e+01,  1.38172700e-02]]],
[[[ -5.35577630e-01,  1.95558600e-02,  1.31108757e+01],
[  1.95558600e-02, -1.51342210e-01, -6.67615000e-03],
[  1.31108757e+01, -6.67615000e-03,  6.90486240e-01]],
[[  2.58256300e-02, -1.49072770e-01, -7.38843000e-03],
[ -1.49072770e-01, -3.63410500e-02,  1.32039847e+01],
[ -7.38843000e-03,  1.32039847e+01,  1.38172700e-02]],
[[  1.33639532e+01, -1.26331100e-02,  6.84650400e-01],
[ -1.26331100e-02,  1.34222177e+01,  1.67851800e-02],
[  6.84650400e-01,  1.67851800e-02,  4.89151396e+01]]]])

g = np.array([[ 0.79389393,  0.54184237,  0.27593346],
[-0.59925749,  0.62028664,  0.50609776],
[ 0.10306737, -0.56714313,  0.8171449 ]])

for i in range(100):
Tprime = rotT(T,g)
``````

``````def rotT(T, g):
Tprime = T
for i in range(4):
slices = [None] * 4
slices[i] = slice(None)
slices *= 2
Tprime = g[slices].T * Tprime
return Tprime.sum(-1).sum(-1).sum(-1).sum(-1)
``````

``````>>> Trot = np.einsum('ai,bj,ck,dl,abcd->ijkl', g, g, g, g, T)
``````

``````import numpy as np
import parakeet

@parakeet.jit
def rotT(T, g):
# ...
``````

``````andrew      10 loops, best of 3: 206 msec per loop
andrew_jit  10 loops, best of 3: 13.3 msec per loop
sven        100 loops, best of 3: 2.39 msec per loop
philipp     1000 loops, best of 3: 0.879 msec per loop
``````