[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 10:54:35 EDT 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