Back to index

python-biopython  1.60
Vector.py
Go to the documentation of this file.
00001 # Copyright (C) 2004, Thomas Hamelryck (thamelry@binf.ku.dk)
00002 # This code is part of the Biopython distribution and governed by its
00003 # license.  Please see the LICENSE file that should have been included
00004 # as part of this package.
00005 
00006 """Vector class, including rotation-related functions."""
00007 
00008 import numpy
00009 
00010 
00011 def m2rotaxis(m):
00012     """
00013     Return angles, axis pair that corresponds to rotation matrix m.
00014     """
00015     # Angle always between 0 and pi
00016     # Sense of rotation is defined by axis orientation
00017     t=0.5*(numpy.trace(m)-1)
00018     t=max(-1, t)
00019     t=min(1, t)
00020     angle=numpy.arccos(t)
00021     if angle<1e-15:
00022         # Angle is 0
00023         return 0.0, Vector(1,0,0)
00024     elif angle<numpy.pi:
00025         # Angle is smaller than pi
00026         x=m[2,1]-m[1,2]
00027         y=m[0,2]-m[2,0]
00028         z=m[1,0]-m[0,1]
00029         axis=Vector(x,y,z)
00030         axis.normalize()
00031         return angle, axis
00032     else:
00033         # Angle is pi - special case!
00034         m00=m[0,0]
00035         m11=m[1,1]
00036         m22=m[2,2]
00037         if m00>m11 and m00>m22:
00038             x=numpy.sqrt(m00-m11-m22+0.5)
00039             y=m[0,1]/(2*x)
00040             z=m[0,2]/(2*x)
00041         elif m11>m00 and m11>m22:
00042             y=numpy.sqrt(m11-m00-m22+0.5)
00043             x=m[0,1]/(2*y)
00044             z=m[1,2]/(2*y)
00045         else:
00046             z=numpy.sqrt(m22-m00-m11+0.5)
00047             x=m[0,2]/(2*z)
00048             y=m[1,2]/(2*z)
00049         axis=Vector(x,y,z)
00050         axis.normalize()
00051         return numpy.pi, axis
00052 
00053 
00054 def vector_to_axis(line, point):
00055     """
00056     Returns the vector between a point and
00057     the closest point on a line (ie. the perpendicular
00058     projection of the point on the line).
00059 
00060     @type line: L{Vector}
00061     @param line: vector defining a line
00062 
00063     @type point: L{Vector}
00064     @param point: vector defining the point
00065     """
00066     line=line.normalized()
00067     np=point.norm()
00068     angle=line.angle(point)
00069     return point-line**(np*numpy.cos(angle))
00070 
00071 
00072 def rotaxis2m(theta, vector):
00073     """
00074     Calculate a left multiplying rotation matrix that rotates
00075     theta rad around vector.
00076 
00077     Example: 
00078     
00079         >>> m=rotaxis(pi, Vector(1,0,0))
00080         >>> rotated_vector=any_vector.left_multiply(m)
00081 
00082     @type theta: float
00083     @param theta: the rotation angle
00084 
00085 
00086     @type vector: L{Vector}
00087     @param vector: the rotation axis
00088 
00089     @return: The rotation matrix, a 3x3 Numeric array.
00090     """
00091     vector=vector.copy()
00092     vector.normalize()
00093     c=numpy.cos(theta)
00094     s=numpy.sin(theta)
00095     t=1-c
00096     x,y,z=vector.get_array()
00097     rot=numpy.zeros((3,3))
00098     # 1st row
00099     rot[0,0]=t*x*x+c
00100     rot[0,1]=t*x*y-s*z
00101     rot[0,2]=t*x*z+s*y
00102     # 2nd row
00103     rot[1,0]=t*x*y+s*z
00104     rot[1,1]=t*y*y+c
00105     rot[1,2]=t*y*z-s*x
00106     # 3rd row
00107     rot[2,0]=t*x*z-s*y
00108     rot[2,1]=t*y*z+s*x
00109     rot[2,2]=t*z*z+c
00110     return rot
00111 
00112 rotaxis=rotaxis2m
00113 
00114 def refmat(p,q):
00115     """
00116     Return a (left multiplying) matrix that mirrors p onto q.
00117 
00118     Example:
00119         >>> mirror=refmat(p,q)
00120         >>> qq=p.left_multiply(mirror)
00121         >>> print q, qq # q and qq should be the same
00122 
00123     @type p,q: L{Vector}
00124     @return: The mirror operation, a 3x3 Numeric array. 
00125     """
00126     p.normalize()
00127     q.normalize()
00128     if (p-q).norm()<1e-5:
00129         return numpy.identity(3)
00130     pq=p-q
00131     pq.normalize()
00132     b=pq.get_array()
00133     b.shape=(3, 1)
00134     i=numpy.identity(3)
00135     ref=i-2*numpy.dot(b, numpy.transpose(b))
00136     return ref
00137 
00138 def rotmat(p,q):
00139     """
00140     Return a (left multiplying) matrix that rotates p onto q.
00141 
00142     Example:
00143         >>> r=rotmat(p,q)
00144         >>> print q, p.left_multiply(r)
00145 
00146     @param p: moving vector
00147     @type p: L{Vector}
00148 
00149     @param q: fixed vector
00150     @type q: L{Vector}
00151 
00152     @return: rotation matrix that rotates p onto q
00153     @rtype: 3x3 Numeric array
00154     """
00155     rot=numpy.dot(refmat(q, -p), refmat(p, -p))
00156     return rot
00157 
00158 def calc_angle(v1, v2, v3):
00159     """
00160     Calculate the angle between 3 vectors
00161     representing 3 connected points.
00162 
00163     @param v1, v2, v3: the tree points that define the angle
00164     @type v1, v2, v3: L{Vector}
00165     
00166     @return: angle
00167     @rtype: float
00168     """
00169     v1=v1-v2
00170     v3=v3-v2
00171     return v1.angle(v3)
00172 
00173 def calc_dihedral(v1, v2, v3, v4):
00174     """
00175     Calculate the dihedral angle between 4 vectors
00176     representing 4 connected points. The angle is in
00177     ]-pi, pi].
00178 
00179     @param v1, v2, v3, v4: the four points that define the dihedral angle
00180     @type v1, v2, v3, v4: L{Vector}
00181     """
00182     ab=v1-v2
00183     cb=v3-v2
00184     db=v4-v3
00185     u=ab**cb
00186     v=db**cb
00187     w=u**v
00188     angle=u.angle(v)
00189     # Determine sign of angle
00190     try:
00191         if cb.angle(w)>0.001:
00192             angle=-angle
00193     except ZeroDivisionError:
00194         # dihedral=pi
00195         pass
00196     return angle
00197 
00198 class Vector(object):
00199     "3D vector"
00200 
00201     def __init__(self, x, y=None, z=None):
00202         if y is None and z is None:
00203             # Array, list, tuple...
00204             if len(x)!=3:
00205                 raise ValueError("Vector: x is not a "
00206                                  "list/tuple/array of 3 numbers")
00207             self._ar=numpy.array(x, 'd')
00208         else:
00209             # Three numbers
00210             self._ar=numpy.array((x, y, z), 'd')
00211 
00212     def __repr__(self):
00213         x,y,z=self._ar
00214         return "<Vector %.2f, %.2f, %.2f>" % (x,y,z)
00215 
00216     def __neg__(self):
00217         "Return Vector(-x, -y, -z)"
00218         a=-self._ar
00219         return Vector(a)
00220 
00221     def __add__(self, other):
00222         "Return Vector+other Vector or scalar"
00223         if isinstance(other, Vector):
00224             a=self._ar+other._ar
00225         else:
00226             a=self._ar+numpy.array(other)
00227         return Vector(a)
00228 
00229     def __sub__(self, other):
00230         "Return Vector-other Vector or scalar"
00231         if isinstance(other, Vector):
00232             a=self._ar-other._ar
00233         else:
00234             a=self._ar-numpy.array(other)
00235         return Vector(a)
00236 
00237     def __mul__(self, other):
00238         "Return Vector.Vector (dot product)"
00239         return sum(self._ar*other._ar)
00240 
00241     def __div__(self, x):
00242         "Return Vector(coords/a)"
00243         a=self._ar/numpy.array(x)
00244         return Vector(a)
00245 
00246     def __pow__(self, other):
00247         "Return VectorxVector (cross product) or Vectorxscalar"
00248         if isinstance(other, Vector):
00249             a,b,c=self._ar
00250             d,e,f=other._ar
00251             c1=numpy.linalg.det(numpy.array(((b,c), (e,f))))
00252             c2=-numpy.linalg.det(numpy.array(((a,c), (d,f))))
00253             c3=numpy.linalg.det(numpy.array(((a,b), (d,e))))
00254             return Vector(c1,c2,c3)
00255         else:
00256             a=self._ar*numpy.array(other)
00257             return Vector(a)
00258 
00259     def __getitem__(self, i):
00260         return self._ar[i]
00261 
00262     def __setitem__(self, i, value):
00263         self._ar[i]=value
00264 
00265     def __contains__(self, i):
00266         return (i in self._ar)
00267 
00268     def norm(self):
00269         "Return vector norm"
00270         return numpy.sqrt(sum(self._ar*self._ar))
00271 
00272     def normsq(self):
00273         "Return square of vector norm"
00274         return abs(sum(self._ar*self._ar))
00275 
00276     def normalize(self):
00277         "Normalize the Vector"
00278         self._ar=self._ar/self.norm()
00279 
00280     def normalized(self):
00281         "Return a normalized copy of the Vector"
00282         v=self.copy()
00283         v.normalize()
00284         return v
00285 
00286     def angle(self, other):
00287         "Return angle between two vectors"
00288         n1=self.norm()
00289         n2=other.norm()
00290         c=(self*other)/(n1*n2)
00291         # Take care of roundoff errors
00292         c=min(c,1)
00293         c=max(-1,c)
00294         return numpy.arccos(c)
00295 
00296     def get_array(self):
00297         "Return (a copy of) the array of coordinates"
00298         return numpy.array(self._ar)
00299 
00300     def left_multiply(self, matrix):
00301         "Return Vector=Matrix x Vector"
00302         a=numpy.dot(matrix, self._ar)
00303         return Vector(a)
00304 
00305     def right_multiply(self, matrix):
00306         "Return Vector=Vector x Matrix"
00307         a=numpy.dot(self._ar, matrix)
00308         return Vector(a)
00309 
00310     def copy(self):
00311         "Return a deep copy of the Vector"
00312         return Vector(self._ar)
00313 
00314 if __name__=="__main__":
00315 
00316         from numpy.random import random
00317 
00318         v1=Vector(0,0,1)
00319         v2=Vector(0,0,0)
00320         v3=Vector(0,1,0)
00321         v4=Vector(1,1,0)
00322 
00323         v4.normalize()
00324 
00325         print v4
00326 
00327         print calc_angle(v1, v2, v3)
00328         dih=calc_dihedral(v1, v2, v3, v4)
00329         # Test dihedral sign
00330         assert(dih>0)
00331         print "DIHEDRAL ", dih
00332 
00333         ref=refmat(v1, v3)
00334         rot=rotmat(v1, v3)
00335 
00336         print v3
00337         print v1.left_multiply(ref)
00338         print v1.left_multiply(rot)
00339         print v1.right_multiply(numpy.transpose(rot))
00340 
00341         # -
00342         print v1-v2
00343         print v1-1
00344         print v1+(1,2,3)
00345         # +
00346         print v1+v2
00347         print v1+3
00348         print v1-(1,2,3)
00349         # *
00350         print v1*v2
00351         # /
00352         print v1/2
00353         print v1/(1,2,3)
00354         # **
00355         print v1**v2
00356         print v1**2
00357         print v1**(1,2,3)
00358         # norm
00359         print v1.norm()
00360         # norm squared
00361         print v1.normsq()
00362         # setitem
00363         v1[2]=10
00364         print v1
00365         # getitem
00366         print v1[2]
00367 
00368         print numpy.array(v1)
00369 
00370         print "ROT"
00371 
00372         angle=random()*numpy.pi
00373         axis=Vector(random(3)-random(3))
00374         axis.normalize()
00375 
00376         m=rotaxis(angle, axis)
00377 
00378         cangle, caxis=m2rotaxis(m)
00379 
00380         print angle-cangle
00381         print axis-caxis
00382         print