Welche Vorschläge haben die Leute, um Diederwinkel in Python schnell zu berechnen?

In den Diagrammen ist phi der Diederwinkel:

Dihedral angle1

Dihedral angle2

Was ist Ihr Bestes für die Berechnung von Winkeln im Bereich von 0 bis pi? Was ist mit 0 bis 2pi?

"Best" bedeutet hier eine Mischung aus schnell und numerisch stabil. Methoden, die Werte über den gesamten Bereich von 0 bis 2 pi zurückgeben, werden bevorzugt, aber wenn Sie eine unglaublich schnelle Methode zur Berechnung der Dieder über 0 haben, teilen Sie dies auch.

Hier sind meine 3 besten Bemühungen. Nur der zweite gibt Winkel zwischen 0 und 2 pi zurück. Es ist auch das langsamste.

Allgemeine Kommentare zu meinen Ansätzen:

Arccos () in Numpy scheint ziemlich stabil zu sein, aber da die Leute dieses Problem ansprechen, verstehe ich es möglicherweise nicht ganz.

Die Verwendung von einsum kam von hier. Warum ist numpys einsum schneller als numpys eingebaute Funktionen?

Die Diagramme und einige Inspirationen kamen von hier. Wie berechne ich einen Diederwinkel bei kartesischen Koordinaten? ?

Die 3 Ansätze mit Kommentaren:

import numpy as np
from time import time

# This approach tries to minimize magnitude and sqrt calculations
def dihedral1(p):
    # Calculate vectors between points, b1, b2, and b3 in the diagram
    b = p[:-1] - p[1:]
    # "Flip" the first vector so that eclipsing vectors have dihedral=0
    b[0] *= -1
    # Use dot product to find the components of b1 and b3 that are not
    # perpendicular to b2. Subtract those components. The resulting vectors
    # lie in parallel planes.
    v = np.array( [ v - (v.dot(b[1])/b[1].dot(b[1])) * b[1] for v in [b[0], b[2]] ] )
    # Use the relationship between cos and dot product to find the desired angle.
    return np.degrees(np.arccos( v[0].dot(v[1])/(np.linalg.norm(v[0]) * np.linalg.norm(v[1]))))

# This is the straightforward approach as outlined in the answers to
# "How do I calculate a dihedral angle given Cartesian coordinates?"
def dihedral2(p):
    b = p[:-1] - p[1:]
    b[0] *= -1
    v = np.array( [ v - (v.dot(b[1])/b[1].dot(b[1])) * b[1] for v in [b[0], b[2]] ] )
    # Normalize vectors
    v /= np.sqrt(np.einsum('...i,...i', v, v)).reshape(-1,1)
    b1 = b[1] / np.linalg.norm(b[1])
    x = np.dot(v[0], v[1])
    m = np.cross(v[0], b1)
    y = np.dot(m, v[1])
    return np.degrees(np.arctan2( y, x ))

# This one starts with two cross products to get a vector perpendicular to
# b2 and b1 and another perpendicular to b2 and b3. The angle between those vectors
# is the dihedral angle.
def dihedral3(p):
    b = p[:-1] - p[1:]
    b[0] *= -1
    v = np.array( [np.cross(v,b[1]) for v in [b[0], b[2]] ] )
    # Normalize vectors
    v /= np.sqrt(np.einsum('...i,...i', v, v)).reshape(-1,1)
    return np.degrees(np.arccos( v[0].dot(v[1]) ))

dihedrals = [ ("dihedral1", dihedral1), ("dihedral2", dihedral2), ("dihedral3", dihedral3) ]

Benchmarking:

# Testing arccos near 0
# Answer is 0.000057
p1 = np.array([
                [ 1,           0,         0     ],
                [ 0,           0,         0     ],
                [ 0,           0,         1     ],
                [ 0.999999,    0.000001,  1     ]
                ])

# +x,+y
p2 = np.array([
                [ 1,           0,         0     ],
                [ 0,           0,         0     ],
                [ 0,           0,         1     ],
                [ 0.1,         0.6,       1     ]
                ])

# -x,+y
p3 = np.array([
                [ 1,           0,         0     ],
                [ 0,           0,         0     ],
                [ 0,           0,         1     ],
                [-0.3,         0.6,       1     ]
                ])
# -x,-y
p4 = np.array([
                [ 1,           0,         0     ],
                [ 0,           0,         0     ],
                [ 0,           0,         1     ],
                [-0.3,        -0.6,       1     ]
                ])
# +x,-y
p5 = np.array([
                [ 1,           0,         0     ],
                [ 0,           0,         0     ],
                [ 0,           0,         1     ],
                [ 0.6,        -0.6,       1     ]
                ])

for d in dihedrals:
    name = d[0]
    f = d[1]
    print "%s:    %12.6f    %12.6f    %12.6f    %12.6f    %12.6f" \
            % (name, f(p1), f(p2), f(p3), f(p4), f(p5))
print

def profileDihedrals(f):
    t0 = time()
    for i in range(20000):
        p = np.random.random( (4,3) )
        f(p)
        p = np.random.randn( 4,3 )
        f(p)
    return(time() - t0)

print "dihedral1:    ", profileDihedrals(dihedral1)
print "dihedral2:    ", profileDihedrals(dihedral2)
print "dihedral3:    ", profileDihedrals(dihedral3)

Benchmarking-Ausgabe:

dihedral1:        0.000057       80.537678      116.565051      116.565051       45.000000
dihedral2:        0.000057       80.537678      116.565051     -116.565051      -45.000000
dihedral3:        0.000057       80.537678      116.565051      116.565051       45.000000

dihedral1:     2.79781794548
dihedral2:     3.74271392822
dihedral3:     2.49604296684

Wie Sie im Benchmarking sehen können, ist der letzte tendenziell der schnellste, während der zweite der einzige ist, der Winkel aus dem gesamten Bereich von 0 bis 2 pi zurückgibt, da er arctan2 verwendet.

28
Praxeolitic 1 Dez. 2013 im 00:28

3 Antworten

Beste Antwort

Hier ist eine Implementierung für den Torsionswinkel über den gesamten 2pi-Bereich, die etwas schneller ist, nicht auf numpy Macken zurückgreift (einsum ist auf mysteriöse Weise schneller als logisch äquivalenter Code) und einfacher zu lesen ist.

Hier gibt es sogar ein bisschen mehr als nur Hacks - auch die Mathematik ist anders. Die in dihedral2 der Frage verwendete Formel verwendet 3 Quadratwurzeln und 1 Kreuzprodukt, die Formel auf Wikipedia verwendet 1 Quadratwurzel und 3 Kreuzprodukte, aber die in der folgenden Funktion verwendete Formel verwendet nur 1 Quadratwurzel und 1 Kreuzprodukt . Dies ist wahrscheinlich so einfach wie die Mathematik nur sein kann.

Funktionen mit 2pi-Bereichsfunktion aus Frage, Wikipedia-Vergleichsformel und der neuen Funktion:

Dihedrals.py

#!/usr/bin/env python
# -*- coding: utf-8 -*-

import numpy as np

def old_dihedral2(p):
    """http://stackoverflow.com/q/20305272/1128289"""
    b = p[:-1] - p[1:]
    b[0] *= -1
    v = np.array( [ v - (v.dot(b[1])/b[1].dot(b[1])) * b[1] for v in [b[0], b[2]] ] )
    # Normalize vectors
    v /= np.sqrt(np.einsum('...i,...i', v, v)).reshape(-1,1)
    b1 = b[1] / np.linalg.norm(b[1])
    x = np.dot(v[0], v[1])
    m = np.cross(v[0], b1)
    y = np.dot(m, v[1])
    return np.degrees(np.arctan2( y, x ))


def wiki_dihedral(p):
    """formula from Wikipedia article on "Dihedral angle"; formula was removed
    from the most recent version of article (no idea why, the article is a
    mess at the moment) but the formula can be found in at this permalink to
    an old version of the article:
    https://en.wikipedia.org/w/index.php?title=Dihedral_angle&oldid=689165217#Angle_between_three_vectors
    uses 1 sqrt, 3 cross products"""
    p0 = p[0]
    p1 = p[1]
    p2 = p[2]
    p3 = p[3]

    b0 = -1.0*(p1 - p0)
    b1 = p2 - p1
    b2 = p3 - p2

    b0xb1 = np.cross(b0, b1)
    b1xb2 = np.cross(b2, b1)

    b0xb1_x_b1xb2 = np.cross(b0xb1, b1xb2)

    y = np.dot(b0xb1_x_b1xb2, b1)*(1.0/np.linalg.norm(b1))
    x = np.dot(b0xb1, b1xb2)

    return np.degrees(np.arctan2(y, x))


def new_dihedral(p):
    """Praxeolitic formula
    1 sqrt, 1 cross product"""
    p0 = p[0]
    p1 = p[1]
    p2 = p[2]
    p3 = p[3]

    b0 = -1.0*(p1 - p0)
    b1 = p2 - p1
    b2 = p3 - p2

    # normalize b1 so that it does not influence magnitude of vector
    # rejections that come next
    b1 /= np.linalg.norm(b1)

    # vector rejections
    # v = projection of b0 onto plane perpendicular to b1
    #   = b0 minus component that aligns with b1
    # w = projection of b2 onto plane perpendicular to b1
    #   = b2 minus component that aligns with b1
    v = b0 - np.dot(b0, b1)*b1
    w = b2 - np.dot(b2, b1)*b1

    # angle between v and w in a plane is the torsion angle
    # v and w may not be normalized but that's fine since tan is y/x
    x = np.dot(v, w)
    y = np.dot(np.cross(b1, v), w)
    return np.degrees(np.arctan2(y, x))

Die neue Funktion würde wahrscheinlich mit 4 separaten Argumenten etwas bequemer aufgerufen, aber um der Signatur in der ursprünglichen Frage zu entsprechen, wird das Argument einfach sofort entpackt.

Code zum Testen:

Test_dihedrals.ph

from dihedrals import *

# some atom coordinates for testing
p0 = np.array([24.969, 13.428, 30.692]) # N
p1 = np.array([24.044, 12.661, 29.808]) # CA
p2 = np.array([22.785, 13.482, 29.543]) # C
p3 = np.array([21.951, 13.670, 30.431]) # O
p4 = np.array([23.672, 11.328, 30.466]) # CB
p5 = np.array([22.881, 10.326, 29.620]) # CG
p6 = np.array([23.691,  9.935, 28.389]) # CD1
p7 = np.array([22.557,  9.096, 30.459]) # CD2

# I guess these tests do leave 1 quadrant (-x, +y) untested, oh well...

def test_old_dihedral2():
    assert(abs(old_dihedral2(np.array([p0, p1, p2, p3])) - (-71.21515)) < 1E-4)
    assert(abs(old_dihedral2(np.array([p0, p1, p4, p5])) - (-171.94319)) < 1E-4)
    assert(abs(old_dihedral2(np.array([p1, p4, p5, p6])) - (60.82226)) < 1E-4)
    assert(abs(old_dihedral2(np.array([p1, p4, p5, p7])) - (-177.63641)) < 1E-4)


def test_new_dihedral1():
    assert(abs(wiki_dihedral(np.array([p0, p1, p2, p3])) - (-71.21515)) < 1E-4)
    assert(abs(wiki_dihedral(np.array([p0, p1, p4, p5])) - (-171.94319)) < 1E-4)
    assert(abs(wiki_dihedral(np.array([p1, p4, p5, p6])) - (60.82226)) < 1E-4)
    assert(abs(wiki_dihedral(np.array([p1, p4, p5, p7])) - (-177.63641)) < 1E-4)


def test_new_dihedral2():
    assert(abs(new_dihedral(np.array([p0, p1, p2, p3])) - (-71.21515)) < 1E-4)
    assert(abs(new_dihedral(np.array([p0, p1, p4, p5])) - (-171.94319)) < 1E-4)
    assert(abs(new_dihedral(np.array([p1, p4, p5, p6])) - (60.82226)) < 1E-4)
    assert(abs(new_dihedral(np.array([p1, p4, p5, p7])) - (-177.63641)) < 1E-4)

Code für das Timing:

Time_dihedrals.py

#!/usr/bin/env python
# -*- coding: utf-8 -*-

from dihedrals import *
from time import time

def profileDihedrals(f):
    t0 = time()
    for i in range(20000):
        p = np.random.random( (4,3) )
        f(p)
        p = np.random.randn( 4,3 )
        f(p)
    return(time() - t0)

print("old_dihedral2: ", profileDihedrals(old_dihedral2))
print("wiki_dihedral: ", profileDihedrals(wiki_dihedral))
print("new_dihedral:  ", profileDihedrals(new_dihedral))

Die Funktionen können mit pytest als pytest ./test_dihedrals.py getestet werden.

Timing-Ergebnisse:

./time_dihedrals.py
old_dihedral2: 1.6442952156066895
wiki_dihedral: 1.3895585536956787
new_dihedral:  0.8703620433807373

new_dihedral ist ungefähr doppelt so schnell wie old_dihedral2.

... Sie können auch sehen, dass die für diese Antwort verwendete Hardware viel leistungsfähiger ist als die in der Frage verwendete Hardware (3.74 vs 1.64 für dihedral2) ;-P

Wenn Sie noch aggressiver werden möchten, können Sie Pypy verwenden. Zum Zeitpunkt des Schreibens unterstützt pypy numpy.cross nicht, aber Sie können stattdessen einfach ein in Python implementiertes Cross-Produkt verwenden. Für ein 3-Vektor-Kreuzprodukt ist der C-Pypy wahrscheinlich mindestens so gut wie das, was Numpy verwendet. Dadurch wird die Zeit für mich auf 0,60 gesenkt, aber dabei waten wir in dummen Hax.

Gleicher Benchmark, aber mit der gleichen Hardware wie in der Frage verwendet:

old_dihedral2: 3.0171279907226562
wiki_dihedral: 3.415065050125122
new_dihedral:  2.086946964263916
17
Praxeolitic 6 Mai 2016 im 19:34

Hier ist die endgültige Antwort. Den Autoren stehen sowohl Python-Versionen als auch die C-Version zur Verfügung.

http://onlinelibrary.wiley.com/doi/10.1002/jcc.20237/abstract

Praktische Umwandlung vom Torsionsraum in den kartesischen Raum für die In-Silico-Proteinsynthese

First published: 16 May 2005
-1
Charlie Strauss 2 Sept. 2016 im 16:56

Mein Ansatz:

Bilden Sie die Vektoren b4 = b1 / \ b2 und b5 = b2 / \ b4. Diese bilden mit b2 einen orthogonalen Rahmen und die Länge von b5 ist die von b2 mal der von b4 (da sie orthogonal sind). Wenn Sie b3 auf diese beiden Vektoren projizieren, erhalten Sie (skalierte) 2D-Koordinaten der Spitze von b3 wie in Ihrer zweiten Figur. Der Winkel wird durch atan2 im Bereich -Pi .. + Pi angegeben.

b4= cross(b1, b2);
b5= cross(b2, b4);
Dihedral= atan2(dot(b3, b4), dot(b3, b5) * sqrt(dot(b2, b2)));

Ähnlich wie bei Ihrer Dieder2. 12 Additionen, 21 Multiplikationen, 1 Quadratwurzel, 1 Arkustangens. Sie können den Ausdruck für b5 mithilfe der Ausweisungsformel umschreiben, dies hilft jedoch nicht wirklich.

VORSICHT: Ich habe die Zeichen- / Quadrantenprobleme nicht gründlich überprüft!

1
Yves Daoust 18 Dez. 2013 im 11:01