functions - numpy python tutorial
Ángulo diédrico/de torsión desde cuatro puntos en coordenadas cartesianas en Python (3)
Aquí está la respuesta definitiva. Los autores tienen versiones de python disponibles, así como la versión C.
http://onlinelibrary.wiley.com/doi/10.1002/jcc.20237/abstract
Conversión práctica del espacio de torsión al espacio cartesiano para la síntesis de proteínas in silico
First published: 16 May 2005
¿Qué sugerencias tienen las personas para calcular rápidamente los ángulos diedros en Python?
En los diagramas, phi es el ángulo diedro:
¿Cuál es tu mejor para calcular ángulos en el rango de 0 a pi? ¿Qué hay de 0 a 2pi?
"Mejor" aquí significa una mezcla de rápida y numéricamente estable. Se prefieren los métodos que devuelven valores en todo el rango de 0 a 2pi, pero si tiene una manera increíblemente rápida de calcular el diedro sobre 0 a pi, también comparta eso.
Aquí están mis 3 mejores esfuerzos. Solo el 2do devuelve ángulos entre 0 y 2pi. También es el más lento.
Comentarios generales sobre mis enfoques:
Arccos () en Numpy parece bastante estable, pero como la gente plantea este problema, es posible que no lo entienda completamente.
El uso de einsum vino de aquí. ¿Por qué ensum de numpy es más rápido que las funciones integradas de numpy?
Los diagramas y algo de inspiración vinieron de aquí. ¿Cómo calculo un ángulo diedro dadas las coordenadas cartesianas?
Los 3 enfoques con comentarios:
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)
Resultados de referencia:
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
Como puede ver en la evaluación comparativa, el último tiende a ser el más rápido, mientras que el segundo es el único que devuelve ángulos desde el rango completo de 0 a 2pi, ya que utiliza arctan2.
Aquí hay una implementación para el ángulo de torsión en el rango completo de 2pi que es un poco más rápido, no recurre a las extravagantes (einsum es misteriosamente más rápido que el código lógicamente equivalente) y es más fácil de leer.
Incluso hay un poco más que simples piratas informáticos, las matemáticas también son diferentes. La fórmula utilizada en dihedral2 la pregunta usa 3 raíces cuadradas y 1 producto cruzado, la fórmula en Wikipedia usa 1 raíz cuadrada y 3 productos cruzados, pero la fórmula utilizada en la siguiente función usa solo 1 raíz cuadrada y 1 producto cruzado. Esto es probablemente tan simple como las matemáticas pueden obtener.
Funciones con rango de 2pi función de pregunta, fórmula de Wikipedia para comparación y la nueva función:
dihedrals.py
#!/usr/bin/env python
# -*- coding: utf-8 -*-
import numpy as np
def old_dihedral2(p):
"""http://.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))
La nueva función probablemente se llame un poco más convenientemente con 4 argumentos separados, pero para que coincida con la firma en la pregunta original, simplemente desempaqueta el argumento inmediatamente.
Código de prueba:
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)
Código de tiempo:
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))
Las funciones se pueden probar con pytest como pytest ./test_dihedrals.py .
Resultados de tiempo:
./time_dihedrals.py
old_dihedral2: 1.6442952156066895
wiki_dihedral: 1.3895585536956787
new_dihedral: 0.8703620433807373
new_dihedral es aproximadamente el doble de rápido que old_dihedral2 .
... también puede ver que el hardware utilizado para esta respuesta es mucho más robusto que el hardware utilizado en la pregunta (3.74 vs 1.64 para dihedral2 ) ;-P
Si quieres ser aún más agresivo puedes usar pypy. En el momento de escribir pypy no es compatible con numpy.cross pero puede usar un producto cruzado implementado en python. Para un producto cruzado de 3 vectores, la generación de C pypy es probablemente al menos tan buena como la que usa numpy. Al hacerlo, el tiempo se reduce a 0,60 para mí, pero en este momento nos estamos metiendo en tonto hax.
El mismo punto de referencia pero con el mismo hardware utilizado en la pregunta:
old_dihedral2: 3.0171279907226562
wiki_dihedral: 3.415065050125122
new_dihedral: 2.086946964263916
Mi acercamiento:
Forme los vectores b4 = b1 / / b2 y b5 = b2 / / b4. Estos forman un marco ortogonal con b2 y la longitud de b5 es la de b2 veces la de b4 (ya que son ortogonales). Proyectar b3 en estos dos vectores le proporciona coordenadas 2D (escaladas) de la punta de b3 como en su segunda figura. El ángulo viene dado por atan2 en el rango -Pi .. + Pi.
b4= cross(b1, b2);
b5= cross(b2, b4);
Dihedral= atan2(dot(b3, b4), dot(b3, b5) * sqrt(dot(b2, b2)));
Similar a tu diedro2. 12 sumas, 21 se multiplica, 1 raíz cuadrada, 1 arctangente. Puedes reescribir la expresión para b5 usando la fórmula de expulsión, pero esto no ayuda.
PRECAUCIÓN: ¡No revisé a fondo los problemas de signos / cuadrantes!