Project Euclidian tensor of fourth order onto subspace of Hooke's tensors

For a given Euclidian tensor of fourth order in three dimensions $mathbf{A}$, how does one compute the Hooke’s part $mathbf{H}$ of $mathbf{A}$?

With the space of Hooke’s tensors $mathcal{H}$, containing tensors with inner symmetry, i.e. both minor symmetries and the major symmetry, i.e. the property $H_{ijkl} = H_{jikl} = H_{ijlk} = H_{klij}$.

[Rychlewski 2000]
introduces the subgroup of permutations
$Sigma_{mathcal{H}} = {(1234), (2134), (1243), (3412), (2143), (4321), (1423), (1324)}$
and a projector
$mathcal{s}_{mathcal{H}} = frac{1}{8}(sum_{i=1}^8 sigma_{i})$ with $sigma_i in Sigma_{mathcal{H}}$
projecting onto the 21-dimensional space of Hooke’s tensors.
This is exactly what is needed.
However, implementing the projector using
$sigma times
(mathbf{a}_1 otimes mathbf{a}_2 otimes mathbf{a}_3 otimes mathbf{a}_4) =
mathbf{a}_{sigma^{-1}(1)} otimes
mathbf{a}_{sigma^{-1}(2)} otimes
mathbf{a}_{sigma^{-1}(3)} otimes
mathbf{a}_{sigma^{-1}(4)}$
with permutation being for example $sigma = (2413)$, meaning $sigma(1) = 2, sigma(2) = 4, sigma(3) = 1, sigma(4) = 3,$ does not yield the expected result.

What have I done wrong or misunderstood?

Implementation:

import numpy as np
A = np.random.rand(3, 3, 3, 3)

permutations_inv = [
    (0, 1, 2, 3),
    (1, 0, 2, 3),
    (0, 1, 3, 2),
    (2, 3, 0, 1),
    (1, 0, 3, 2),
    (3, 2, 1, 0),
    (0, 3, 1, 2),
    (0, 2, 1, 3),
    ]

# Invert permutation according definition
permutations = []
for p in permutations_inv:
    permutations.append([
            p.index(0),
            p.index(1),
            p.index(2),
            p.index(3),
            ])

H = 1.0/len(permutations) * sum(
                            A.transpose(perm)
                            for perm in permutations
                            )


def has_sym_inner(A):
    sym = True
    for i in range(3):
        for j in range(3):
            for k in range(3):
                for l in range(3):
                    if not all(
                                [
                                    A[i, j, k, l] == A[j, i, k, l],
                                    A[i, j, k, l] == A[i, j, l, k],
                                    A[i, j, k, l] == A[k, l, i, j]
                                ]
                                ):
                        sym = False

    return sym

print('has_sym_inner =', has_sym_inner(H))

Verifying visually using
mechkit

import mechkit
con = mechkit.notation.Converter()
print('A=n', con.to_mandel9(A))
print('H=n', con.to_mandel9(H))

[Rychlewski 2000]: Rychlewski, J. (2000). A qualitative approach to Hooke’s tensors. Part I. Archives of Mechanics, 52(4-5), 737-759