Какие предложения у людей есть для быстрого вычисления двугранных углов в Python?
На диаграммах phi - двугранный угол:
Что лучше всего для расчета углов в диапазоне от 0 до pi? Как насчет от 0 до 2pi?
"Лучшее" здесь означает некоторое сочетание быстрого и численно стабильного. Методы, возвращающие значения во всем диапазоне от 0 до 2pi, являются предпочтительными, но если у вас есть невероятно быстрый способ вычисления диэдра в пределах от 0 до pi, то это тоже.
Вот мои 3 лучших усилия. Только второй возвращает углы между 0 и 2pi. Это также самый медленный.
Общие комментарии о моих подходах:
arccos() в Numpy кажется довольно стабильным, но поскольку люди поднимают эту проблему, я могу просто не полностью ее понять.
Использование einsum исходило отсюда. Почему numpy einsum быстрее, чем встроенные функции numpy?
Диаграммы и некоторое вдохновение пришли сюда. Как рассчитать двугранный угол, заданный декартовыми координатами?
3 подхода с комментариями:
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) ]
Бенчмаркинг:
# 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)
Результат тестирования:
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
Как вы можете видеть в бенчмаркинге, последний имеет тенденцию быть самым быстрым, а второй - единственным, который возвращает углы из полного диапазона от 0 до 2pi, поскольку он использует arctan2.
