Эффективное сжатие тензора Леви-Чивита с Numpy эинсумом - PullRequest
Новая
       76

Эффективное сжатие тензора Леви-Чивита с Numpy эинсумом

0 голосов
/ 18 апреля 2020

Я хочу сжать большие n-мерные векторы с тензором Леви Чивиты. Если я хочу использовать функцию Numpy einsum, мне нужно заранее определить тензор Levi Civita, который быстро подорвет память моего компьютера в n измерениях. Мне просто нужно заставить функцию einsum принимать callable, чтобы я мог сделать что-то вроде этого: 

import operator
from functools import reduce

import numpy as np
from math import factorial

def prod(a):
    """
    Return product of elements of a. Start with int 1 so if only ints are included then an int result is
    returned.
    """
    return reduce(operator.mul, a, 1)


def eijk(args):
    """
    Represent the Levi-Civita symbol. For even permutations of indices it returns 1, for odd
    permutations -1, and for everything else (a repeated index) it returns 0. Thus it represents an
    alternating pseudotensor.

    Parameters
    ----------
    args : tuple with int
        A tuple with indices.

    Returns
    -------
    int
    """
    n = len(args)

    return prod(prod(args[j] - args[i] for j in range(i + 1, n)) / factorial(i) for i in range(n))


a = np.array([1., 0., 0., 8.])
b = np.array([1., 5., 8., 4.])

c = np.einsum("ijkl, j, k -> ij", eijk, a, b)

Это приведет к ValueError: einstein sum subscripts string contains too many subscripts for operand 0. Он работает только тогда, когда тензор Levi Civita предопределен следующим образом: 

# The array `eijk_array` has a shape of (4, 4, 4, 4). That's why I won't type the entire tensor. 
eijk_array = array([[[[ 0.,  0.,  0.,  0.],
                      [ 0.,  0.,  0.,  0.],
                      [ 0.,  0.,  0.,  0.],
                      [ 0.,  0.,  0.,  0.]],
                     [[ 0.,  0.,  0.,  0.],
                      [ 0.,  0.,  0.,  0.],
                      [ 0.,  0.,  0.,  1.],
                      [ 0.,  0., -1.,  0.]],
                     [[ 0.,  0.,  0.,  0.],
                      [ 0.,  0.,  0., -1.],
                      [ 0.,  0.,  0.,  0.],
                      [ 0.,  1.,  0.,  0.]],
                     [[ 0.,  0.,  0.,  0.],
                      [ 0.,  0.,  1.,  0.],
                      [ 0., -1.,  0.,  0.],
                      [ 0.,  0.,  0.,  0.]]],
                      ...

a = np.array([1., 0., 0., 8.])
b = np.array([1., 5., 8., 4.])

c = np.einsum("ijkl, j, k -> ij", eijk_array, a, b)

Работа с четырьмя измерениями все еще идет хорошо. Но при работе с 10, 20 или 30 размерными векторами, как в моем случае, он не может быть рассчитан в этой форме. Для трехмерного вектора я решил это, как показано в следующем примере: 

a = np.array([1., 0., 0.])
b = np.array([1., 5., 8.])

args = [a, b]

# The dimension of the output. This is the amount of `.. -> ij` `np.einsum` subscriptions 
dim = a.shape[0] - len(args)

# Start Contraction (`eijk` is the predefined function above)
c = np.zeros(np.repeat(a.shape[0], dim))

for k in range(a.shape[0]):
    for i, j in product(range(a.shape[0]), repeat=2):
        c[k] = c[k] + a[i] * b[j] * eijk((i, j, k))

Я пока не знаю, как расширить это l oop до n измерений. Прежде чем пытаться решить проблему с большими усилиями, я хотел спросить, есть ли уже решение для этой проблемы.

Редактировать

Я пытался обмануть Numpy, создав подкласс Numpy .ndarray и изменив метод magi c __getitem__, например: 

class LeviCivita(np.ndarray):
    def __new__(...):
        # ....

    def __getitem__(self, item):
        return eijk(item)

К сожалению, это тоже не сработало.

...