[Biopython-dev] [Bug 2592] New: numpy migration for Bio.PDB.Vector
bugzilla-daemon at portal.open-bio.org
bugzilla-daemon at portal.open-bio.org
Mon Sep 22 14:54:35 UTC 2008
http://bugzilla.open-bio.org/show_bug.cgi?id=2592
Summary: numpy migration for Bio.PDB.Vector
Product: Biopython
Version: Not Applicable
Platform: All
OS/Version: All
Status: NEW
Severity: enhancement
Priority: P2
Component: Other
AssignedTo: biopython-dev at biopython.org
ReportedBy: meesters at uni-mainz.de
see http://lists.open-bio.org/pipermail/biopython/2008-September/004505.html
The code is pretty similar to the original one. I don't mind, if it won't be
used.
Vector.py:
from scipy.linalg import det #determinant
from numpy import allclose, arccos, array, cos, dot, eye, float32, matrix, sin,
\
sqrt, sum, trace, transpose, zeros
from math import acos
def m2rotaxis(m):
"""
Return angles, axis pair that corresponds to rotation matrix m.
"""
# Angle always between 0 and pi
# Sense of rotation is defined by axis orientation
t=0.5*(trace(m)-1)
t=max(-1, t)
t=min(1, t)
angle=acos(t)
if angle<1e-15:
# Angle is 0
return 0.0, Vector(1,0,0)
elif angle<pi:
# Angle is smaller than pi
x=m[2,1]-m[1,2]
y=m[0,2]-m[2,0]
z=m[1,0]-m[0,1]
axis=Vector(x,y,z)
axis.normalize()
return angle, axis
else:
# Angle is pi - special case!
m00=m[0,0]
m11=m[1,1]
m22=m[2,2]
if m00>m11 and m00>m22:
x=sqrt(m00-m11-m22+0.5)
y=m[0,1]/(2*x)
z=m[0,2]/(2*x)
elif m11>m00 and m11>m22:
y=sqrt(m11-m00-m22+0.5)
x=m[0,1]/(2*y)
z=m[1,2]/(2*y)
else:
z=sqrt(m22-m00-m11+0.5)
x=m[0,2]/(2*z)
y=m[1,2]/(2*z)
axis=Vector(x,y,z)
axis.normalize()
return pi, axis
def vector_to_axis(line, point):
"""
Returns the vector between a point and
the closest point on a line (ie. the perpendicular
projection of the point on the line).
@type line: L{Vector}
@param line: vector defining a line
@type point: L{Vector}
@param point: vector defining the point
"""
line=line.normalized()
np=point.norm()
angle=line.angle(point)
return point-line**(np*cos(angle))
def calc_angle(v1, v2, v3):
"""
Calculate the angle between 3 vectors
representing 3 connected points.
@param v1, v2, v3: the tree points that define the angle
@type v1, v2, v3: L{Vector}
@return: angle
@rtype: float
"""
v1=v1-v2
v3=v3-v2
return v1.angle(v3)
def calc_dihedral(v1, v2, v3, v4):
"""
Calculate the dihedral angle between 4 vectors
representing 4 connected points. The angle is in
]-pi, pi].
@param v1, v2, v3, v4: the four points that define the dihedral angle
@type v1, v2, v3, v4: L{Vector}
"""
ab=v1-v2
cb=v3-v2
db=v4-v3
u=ab**cb
v=db**cb
w=u**v
angle=u.angle(v)
# Determine sign of angle
try:
if cb.angle(w)>0.001:
angle=-angle
except ZeroDivisionError:
# dihedral=pi
pass
return angle
def rotaxis(theta, vector):
"""
Calculate a left multiplying rotation matrix that rotates
theta rad around vector.
Example:
>>> m=rotaxis(pi, Vector(1,0,0))
>>> rotated_vector=any_vector.left_multiply(m)
@type theta: float
@param theta: the rotation angle
@type vector: L{Vector}
@param vector: the rotation axis
@return: The rotation matrix, a 3x3 Numeric array.
"""
vector=vector.copy()
vector.normalize()
c=cos(theta)
s=sin(theta)
t=1-c
x,y,z=vector.get_array()
rot=zeros((3,3), "d")
# 1st row
rot[0,0]=t*x*x+c
rot[0,1]=t*x*y-s*z
rot[0,2]=t*x*z+s*y
# 2nd row
rot[1,0]=t*x*y+s*z
rot[1,1]=t*y*y+c
rot[1,2]=t*y*z-s*x
# 3rd row
rot[2,0]=t*x*z-s*y
rot[2,1]=t*y*z+s*x
rot[2,2]=t*z*z+c
return rot
def refmat(p,q):
"""
Return a (left multiplying) matrix that mirrors p onto q.
Example:
>>> mirror=refmat(p,q)
>>> qq=p.left_multiply(mirror)
>>> print q, qq # q and qq should be the same
@type p,q: L{Vector}
@return: The mirror operation, a 3x3 Numeric array.
"""
p.normalize()
q.normalize()
if (p-q).norm()<1e-5:
return eye(3)
pq=p-q
pq.normalize()
b=pq.get_array()
b.shape=(3, 1)
i=eye(3)
ref=i-2* dot(b, transpose(b))
return ref
def rotmat(p,q):
"""
Return a (left multiplying) matrix that rotates p onto q.
Example:
>>> r=rotmat(p,q)
>>> print q, p.left_multiply(r)
@param p: moving vector
@type p: L{Vector}
@param q: fixed vector
@type q: L{Vector}
@return: rotation matrix that rotates p onto q
@rtype: 3x3 Numeric array
"""
rot=refmat(q, -p) * refmat(p, -p).transpose()
return rot
class Vector(object):
"3D vector"
def __init__(self, x, y=None, z=None):
if y is None and z is None:
# Array, list, tuple...
if len(x)!=3:
raise "Vector: x is not a list/tuple/array of 3 numbers"
self._ar=array(x)
else:
# Three numbers
self._ar=array([x, y, z])
def __eq__(self, other):
return allclose(self._ar, other._ar, 0.01)
def __ne__(self, other):
return not self.__eq__(other)
def __repr__(self):
x, y, z = self._ar
return "<Vector %.2f, %.2f, %.2f>" % (x, y, z)
def __neg__(self):
"Return Vector(-x, -y, -z)"
return Vector(-self._ar)
def __add__(self, other):
"Return Vector+other Vector or scalar"
if isinstance(other, Vector):
a=self._ar+other._ar
else:
a=self._ar+array(other)
return Vector(a)
def __sub__(self, other):
"Return Vector-other Vector or scalar"
if isinstance(other, Vector):
a=self._ar-other._ar
else:
a=self._ar-array(other)
return Vector(a)
def __mul__(self, other):
"Return Vector.Vector (dot product)"
return sum(self._ar*other._ar)
def __div__(self, x):
"Return Vector(coords/a)"
a=self._ar/array(x)
return Vector(a)
def __pow__(self, other):
"Return VectorxVector (cross product) or Vectorxscalar"
if isinstance(other, Vector):
a,b,c=self._ar
d,e,f=other._ar
c1=det(array(((b,c), (e,f))))
c2=-det(array(((a,c), (d,f))))
c3=det(array(((a,b), (d,e))))
return Vector(c1,c2,c3)
else:
a=self._ar*array(other)
return Vector(a)
def __getitem__(self, i):
return self._ar[i]
def __setitem__(self, i, value):
self._ar[i]=value
def norm(self):
"Return vector norm"
return sqrt(sum(self._ar*self._ar))
def normsq(self):
"Return square of vector norm"
return abs(sum(self._ar*self._ar))
def normalize(self):
"Normalize the Vector"
self._ar=self._ar/self.norm()
def normalized(self):
"Return a normalized copy of the Vector"
v = self.copy()
v.normalize()
return v
def angle(self, other):
"Return angle between two vectors"
n1=self.norm()
n2=other.norm()
c=(self*other)/(n1*n2)
# Take care of roundoff errors
c=min(c,1)
c=max(-1,c)
return arccos(c)
def get_array(self):
"Return (a copy of) the array of coordinates"
return array(self._ar)
def left_multiply(self, matrix):
"Return Vector=Matrix x Vector"
return Vector(dot(matrix, self._ar))
def right_multiply(self, matrix):
"Return Vector=Vector x Matrix"
return Vector(dot(self._ar, matrix))
def copy(self):
"Return a deep copy of the Vector"
return Vector(self._ar)
#xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx
test_Vector.py:
import unittest
from math import pi, degrees
from numpy import array, allclose, transpose
from Vector import Vector, calc_angle, calc_dihedral, refmat, rotmat, \
rotaxis, vector_to_axis
class TestVectorFunctions(unittest.TestCase):
"""
Vector-class test functions
"""
def setUp(self):
self.v1 = Vector(0,0,1)
self.v2 = Vector(0,0,0)
self.v3 = Vector(0,1,0)
self.v4 = Vector(1,1,0)
self.v5 = Vector(-1,-1,0)
self.ref_array = array([[ 1.0, 0.0, 0.0],
[ 0.0, 0.0, 1.0],
[ 0.0, 1.0, 0.0]])
self.tolerance = 0.001
def test__eq__(self):
"""Vector.__eq__ should return Boolean for equality check on two
Vectors: testing True"""
self.assert_(self.v1 == self.v1)
def test__eq__2(self):
"""Vector.__eq__ should return Boolean for equality check on two
Vectors: testing False"""
self.failIf(self.v1 == self.v2)
def test__ne__(self):
"""Vector.__ne__ should return Boolean for non-equal Vectors: testing
True"""
self.assert_(self.v1 != self.v2)
def test__ne__2(self):
"""Vector.__ne__ should return Boolean for non-equal Vectors: testing
False"""
self.failIf(self.v1 != self.v1)
def test__repr__(self):
"""Vector.__repr__ should return a Vector-object as a nice string"""
self.assertEqual(repr(self.v1), "<Vector 0.00, 0.00, 1.00>")
def test__neg__(self):
"""Vector.__neg__ should Vector-object * -1"""
v = Vector(0,0,-1)
self.assertEqual(-self.v1, v)
def test__add__(self):
"""testing Vector.___add__: Vector + Vector"""
v = Vector(1,1,1)
v2 = Vector(1,1,2)
self.assertEqual(self.v1+v, v2)
def test__add__2(self):
"""testing Vector.___add__: Vector + scalar"""
v = Vector(3,3,4)
self.assertEqual(self.v1+3, v)
def test__add__3(self):
"""testing Vector.___add__: Vector + scalars"""
v = Vector(1,2,4)
self.assertEqual(self.v1+(1,2,3), v)
def test__sub__(self):
"""testing Vector.__sub__(): Vector - Vector"""
self.assertEqual(self.v1-self.v1, self.v2)
def test__sub__2(self):
"""testing Vector.__sub__(): Vector-scalar"""
self.assertEqual(self.v1-1, self.v5)
def test__sub__3(self):
"""testing Vector.__sub__(): Vector-scalars"""
v = Vector(-1,-2,-2)
self.assertEqual(self.v1-(1,2,3), v)
def test__mul__(self):
"""testing Vector.__mul__()"""
self.assertEqual(self.v1 * self.v2, 0)
def test__pow__(self):
"""testing Vector.__pow__()"""
self.assertEqual(self.v1** self.v2, self.v2)
def test__getitem__(self):
"""testing Vector.__getitem__"""
self.assertEqual(self.v1[0], 0)
def test__setitem__(self):
"""testing Vector.__setitem__"""
v = self.v3
v[0] = 1
self.assertEqual(v, self.v4)
def testNorm(self):
"""testing Vector.norm()"""
self.assertEqual(self.v4, self.v4)
def testNormsq(self):
"""testing Vector.normsq()"""
self.assertEqual(self.v4, self.v4)
def testNomalize(self):
"""testing Vector.normalize()"""
self.v4.normalize()
v = Vector(0.71, 0.71, 0.00)
self.assertEqual(self.v4, v)
def testNomalized(self):
"""testing Vector.normalized()"""
self.v4.normalize()
v = Vector(0.71, 0.71, 0.00)
self.assertEqual(self.v4, v)
def testAngle(self):
"""testing Vector.angle()"""
self.assertEqual(degrees(self.v2.angle(self.v1)), 180)
def testGetarray(self):
"""testing Vector.get_array()"""
self.assert_(all(self.v1.get_array() == array((0,0,1))))
def testCopy(self):
"""testing Vector.copy()"""
self.assertEqual(self.v1.copy(), self.v1)
def testCalcangle(self):
"""testing calc_angle()"""
self.assertEqual(degrees(calc_angle(self.v1, self.v2, self.v3)), 90.0)
def testRefmat(self):
"""testing refmat()"""
self.assert_(allclose(refmat(self.v1, self.v3), self.ref_array,
self.tolerance))
def testRotmat(self):
"""testing rotmat()"""
self.assert_(allclose(refmat(self.v1, self.v3), self.ref_array,
self.tolerance))
def testLeftmultiply(self):
"""testing Vector.leftmultiply()"""
self.assertEqual(self.v1.left_multiply(self.ref_array), self.v3)
def testRightmultiply(self):
"""testing Vector.rightmultiply()"""
self.assertEqual(self.v1.right_multiply(transpose(self.ref_array)),
self.v3)
def testRotaxis(self):
"""testing rotaxis()"""
a = array([[ -1.0, 0, 0.0],
[0.0, -1.0, 0.0],
[0.0, 0.0, 1.0]])
self.assert_(allclose(rotaxis(pi, self.v1), a, self.tolerance))
def testVector_to_axis(self):
"""testing vector_to_axis"""
self.assertEqual(vector_to_axis(self.v5, self.v1), self.v1)
def testCalc_dihedral(self):
"""testing calc_dihedral"""
self.assertEqual(degrees(calc_dihedral(self.v1, self.v2, self.v3,
self.v4)), 90)
--
Configure bugmail: http://bugzilla.open-bio.org/userprefs.cgi?tab=email
------- You are receiving this mail because: -------
You are the assignee for the bug, or are watching the assignee.
More information about the Biopython-dev
mailing list