Вычислить знаковую площадь множества треугольников

Я думаю, что проблема не в вычислениях, а в том, как данные организованы в памяти. Если вам не нужно сохранять исходные данные, вы можете попытаться минимизировать количество временных объектов в памяти с помощью функции inplace:

def mymix(p):

    p[0][1] -= p[0][0] # x1 = x1 - x0
    p[0][2] -= p[0][0] # x2 = x2 - x0
    p[1][1] -= p[1][0] # y1 = y1 - y0
    p[1][2] -= p[1][0] # y2 = y2 - y0

    p[0][1] *= p[1][2] # x1 = (x1-x0) * (y2-y0)
    p[0][2] *= p[1][1] # x2 = (x2-x0) * (y1-y0)

    p[0][1] -= p[0][2] # r = (x1-x0) * (y2-y0) - (x2-x0) * (y1-y0)
    p[0][1] /= 2

    return p[0][1]

Вы можете использовать Numba или Cython для оптимизации проблем этого типа. Возможно, лучший способ - вычислить площадь со знаком треугольников напрямую, когда вы получите угловые точки из этих треугольников. Даже в однопоточной версии ваша пропускная способность памяти может быть ограничивающим фактором с этими небольшими прибавлениями и умножениями.

Код

import numba as nb
import numpy 

@nb.njit(fastmath=True)
def mix_nb(p):
    assert p.shape[0]==2
    assert p.shape[1]==3

    res=numpy.empty(p.shape[2],dtype=p.dtype)

    for i in range(p.shape[2]):
      res[i]=(+p[0,2,i] * (p[1,0,i] - p[1,1,i])+ p[0,0,i] * (p[1,1,i] - p[1,2,i])+ p[0,1,i] * (p[1,2,i] - p[1,0,i])) /2.

    return res

#Compile
p= np.random.rand(2, 3, 10000)
A=mix_nb(p)

import perfplot

#Benchmark
perfplot.show(
    setup=lambda n: np.random.rand(2, 3, n),
    kernels=[six, mix, mix2, cross, einsum, einsum2, einsum3,mix_nb],
    n_range=[2**k for k in range(19)],
    logx=True,
    logy=True,
    )

Полученные результаты

Это простая версия Cython, которую я включаю сюда только для полноты картины:

import numpy as np
from libc.stdlib cimport malloc, calloc, realloc, free

def calculate_signed_areas(double[:,:,::1] triangles):
    cdef:
        int i, n
        double area
        double[::1] areas
        double x0, x1, x2
        double y0, y1, y2

    n = triangles.shape[2]
    areas = <double[:n]>malloc(n * sizeof(double))
    for i in range(n):
        x0 = triangles[0,0,i]
        y0 = triangles[1,0,i]
        x1 = triangles[0,1,i]
        y1 = triangles[1,1,i]
        x2 = triangles[0,2,i]
        y2 = triangles[1,2,i]
        area = (
            + x2 * (y0 - y1)
            + x0 * (y1 - y2)
            + x1 * (y2 - y0)
        )/2.0
        areas[i] = area
    return np.asarray(areas)

Я частично согласен с ответом @ TrapII: ​​это проблема использования памяти. Однако я считаю более важным то, как данные хранятся в памяти. Проведем следующий эксперимент:

In [1]: import numpy as np

In [2]: p = np.random.random((1000,3,2))

In [3]: p2 = np.array(np.moveaxis(p, (1, 0), (0, 2)).reshape((6, 1000)), order='C')

In [4]: def mix(p):
   ...:     return (
   ...:           p[:,2,0] * (p[:,0,1] - p[:,1,1])
   ...:         + p[:,0,0] * (p[:,1,1] - p[:,2,1])
   ...:         + p[:,1,0] * (p[:,2,1] - p[:,0,1])
   ...:         ) / 2.0
   ...:     

In [5]: def cross2(p2):
   ...:     e0x, e0y, e1x, e1y, e2x, e2y = p2
   ...:     return 0.5 * ((e1x - e0x)*(e2y - e0y) - (e1y - e0y)*(e2x - e0x))
   ...: 

In [6]: %timeit mix(p)
18.5 µs ± 351 ns per loop (mean ± std. dev. of 7 runs, 100000 loops each)

In [7]: %timeit cross2(p2)
9.03 µs ± 90.9 ns per loop (mean ± std. dev. of 7 runs, 100000 loops each)