Kann NumPy einsum() führt eine Kreuzprodukt zwischen den Segmenten einer Trajektorie

Question

I das Kreuzprodukt von aneinander angrenzenden Segmenten einer Flugbahn durchzuführen (xy-Koordinaten) mit dem folgende Skript:

In [129]: 
def func1(xy, s): 
    size = xy.shape[0]-2*s 
    out = np.zeros(size) 
    for i in range(size): 
     p1, p2 = xy[i], xy[i+s]  #segment 1 
     p3, p4 = xy[i+s], xy[i+2*s] #segment 2 
     out[i] = np.cross(p1-p2, p4-p3) 
    return out 

def func2(xy, s): 
    size = xy.shape[0]-2*s 
    p1 = xy[0:size] 
    p2 = xy[s:size+s] 
    p3 = p2 
    p4 = xy[2*s:size+2*s] 

    tmp1 = p1-p2 
    tmp2 = p4-p3 
    return tmp1[:, 0] * tmp2[:, 1] - tmp2[:, 0] * tmp1[:, 1] 


In [136]: 
xy = np.array([[1,2],[2,3],[3,4],[5,6],[7,8],[2,4],[5,2],[9,9],[1,1]]) 
func2(xy, 2) 

Out[136]: 
array([ 0, -3, 16, 1, 22]) 

func1 ist besonders langsam, weil von der inneren Schleife, so schrieb ich das Cross-Produkt selbst um (func2), das um Grössenordnungen schneller ist.

Ist es möglich, die numpy-einsum-Funktion zu verwenden, um die gleiche Berechnung durchzuführen?

Answer 1

einsum Summen von Produkten berechnet nur, aber man konnte das Kreuzprodukt in eine Summe von Produkten durch Umkehren der Spalten von tmp2 und Änderung des Vorzeichens der ersten Spalte Schuhanzieher:

def func3(xy, s): 
    size = xy.shape[0]-2*s 
    tmp1 = xy[0:size] - xy[s:size+s] 
    tmp2 = xy[2*s:size+2*s] - xy[s:size+s] 
    tmp2 = tmp2[:, ::-1] 
    tmp2[:, 0] *= -1 
    return np.einsum('ij,ij->i', tmp1, tmp2) 

---

Aber func3 ist langsamer als func2.

In [80]: xy = np.tile(xy, (1000, 1)) 

In [104]: %timeit func1(xy, 2) 
10 loops, best of 3: 67.5 ms per loop 

In [105]: %timeit func2(xy, 2) 
10000 loops, best of 3: 73.2 µs per loop 

In [106]: %timeit func3(xy, 2) 
10000 loops, best of 3: 108 µs per loop 

Plausibiltätsprüfung:

In [86]: np.allclose(func1(xy, 2), func3(xy, 2)) 
Out[86]: True 

Ich denke, der Grund, warum func2einsum hier ist, weil die Kosten für die Einstellung der Schleife in einsum für nur 2 Iterationen zu teuer im Vergleich zu nur manuell schlagen schreiben aus der Summe, und die Umkehrung und Multiplikation essen auch einige Zeit.

Answer 2

np.cross ist ein schlaues kleines Biest, das ohne Probleme mit Rundfunk umgehen kann. So können Sie umschreiben Ihre func2 als:

def func2(xy, s): 
    size = xy.shape[0]-2*s 
    p1 = xy[0:size] 
    p2 = xy[s:size+s] 
    p3 = p2 
    p4 = xy[2*s:size+2*s] 
    return np.cross(p1-p2, p4-p3) 

und es wird das richtige Ergebnis produzieren:

>>> func2(xy, 2) 
array([ 0, -3, 16, 1, 22]) 

In der neuesten numpy wird es wahrscheinlich ein bisschen schneller als Ihr Code ausführen, wie es geschrieben wurde minimieren Sie die Erstellung von Zwischen-Arrays. Sie können den Quellcode (reines Python) here betrachten.