Source Code for Module Bio.PDB.Vector'

  1  # Copyright (C) 2004, Thomas Hamelryck (thamelry@binf.ku.dk) 
  2  # This code is part of the Biopython distribution and governed by its 
  3  # license.  Please see the LICENSE file that should have been included 
  4  # as part of this package. 
  5   
  6  """Vector class, including rotation-related functions.""" 
  7   
  8  import numpy 
  9   
 10   
 11 -def m2rotaxis(m): 
 12      """ 
 13      Return angles, axis pair that corresponds to rotation matrix m. 
 14      """ 
 15      # Angle always between 0 and pi 
 16      # Sense of rotation is defined by axis orientation 
 17      t=0.5*(numpy.trace(m)-1) 
 18      t=max(-1, t) 
 19      t=min(1, t) 
 20      angle=numpy.arccos(t) 
 21      if angle<1e-15: 
 22          # Angle is 0 
 23          return 0.0, Vector(1,0,0) 
 24      elif angle<numpy.pi: 
 25          # Angle is smaller than pi 
 26          x=m[2,1]-m[1,2] 
 27          y=m[0,2]-m[2,0] 
 28          z=m[1,0]-m[0,1] 
 29          axis=Vector(x,y,z) 
 30          axis.normalize() 
 31          return angle, axis 
 32      else: 
 33          # Angle is pi - special case! 
 34          m00=m[0,0] 
 35          m11=m[1,1] 
 36          m22=m[2,2] 
 37          if m00>m11 and m00>m22: 
 38              x=numpy.sqrt(m00-m11-m22+0.5) 
 39              y=m[0,1]/(2*x) 
 40              z=m[0,2]/(2*x) 
 41          elif m11>m00 and m11>m22: 
 42              y=numpy.sqrt(m11-m00-m22+0.5) 
 43              x=m[0,1]/(2*y) 
 44              z=m[1,2]/(2*y) 
 45          else: 
 46              z=numpy.sqrt(m22-m00-m11+0.5) 
 47              x=m[0,2]/(2*z) 
 48              y=m[1,2]/(2*z) 
 49          axis=Vector(x,y,z) 
 50          axis.normalize() 
 51          return numpy.pi, axis 
 52   
 53   
 54 -def vector_to_axis(line, point): 
 55      """ 
 56      Returns the vector between a point and 
 57      the closest point on a line (ie. the perpendicular 
 58      projection of the point on the line). 
 59   
 60      @type line: L{Vector} 
 61      @param line: vector defining a line 
 62   
 63      @type point: L{Vector} 
 64      @param point: vector defining the point 
 65      """ 
 66      line=line.normalized() 
 67      np=point.norm() 
 68      angle=line.angle(point) 
 69      return point-line**(np*numpy.cos(angle)) 
 70   
 71   
 72 -def rotaxis2m(theta, vector): 
 73      """ 
 74      Calculate a left multiplying rotation matrix that rotates 
 75      theta rad around vector. 
 76   
 77      Example: 
 78   
 79          >>> m=rotaxis(pi, Vector(1,0,0)) 
 80          >>> rotated_vector=any_vector.left_multiply(m) 
 81   
 82      @type theta: float 
 83      @param theta: the rotation angle 
 84   
 85   
 86      @type vector: L{Vector} 
 87      @param vector: the rotation axis 
 88   
 89      @return: The rotation matrix, a 3x3 Numeric array. 
 90      """ 
 91      vector=vector.copy() 
 92      vector.normalize() 
 93      c=numpy.cos(theta) 
 94      s=numpy.sin(theta) 
 95      t=1-c 
 96      x,y,z=vector.get_array() 
 97      rot=numpy.zeros((3,3)) 
 98      # 1st row 
 99      rot[0,0]=t*x*x+c 
100      rot[0,1]=t*x*y-s*z 
101      rot[0,2]=t*x*z+s*y 
102      # 2nd row 
103      rot[1,0]=t*x*y+s*z 
104      rot[1,1]=t*y*y+c 
105      rot[1,2]=t*y*z-s*x 
106      # 3rd row 
107      rot[2,0]=t*x*z-s*y 
108      rot[2,1]=t*y*z+s*x 
109      rot[2,2]=t*z*z+c 
110      return rot 
111   
112  rotaxis=rotaxis2m 
113   
114   
115 -def refmat(p,q): 
116      """ 
117      Return a (left multiplying) matrix that mirrors p onto q. 
118   
119      Example: 
120          >>> mirror=refmat(p,q) 
121          >>> qq=p.left_multiply(mirror) 
122          >>> print q, qq # q and qq should be the same 
123   
124      @type p,q: L{Vector} 
125      @return: The mirror operation, a 3x3 Numeric array. 
126      """ 
127      p.normalize() 
128      q.normalize() 
129      if (p-q).norm()<1e-5: 
130          return numpy.identity(3) 
131      pq=p-q 
132      pq.normalize() 
133      b=pq.get_array() 
134      b.shape=(3, 1) 
135      i=numpy.identity(3) 
136      ref=i-2*numpy.dot(b, numpy.transpose(b)) 
137      return ref 
138   
139   
140 -def rotmat(p,q): 
141      """ 
142      Return a (left multiplying) matrix that rotates p onto q. 
143   
144      Example: 
145          >>> r=rotmat(p,q) 
146          >>> print q, p.left_multiply(r) 
147   
148      @param p: moving vector 
149      @type p: L{Vector} 
150   
151      @param q: fixed vector 
152      @type q: L{Vector} 
153   
154      @return: rotation matrix that rotates p onto q 
155      @rtype: 3x3 Numeric array 
156      """ 
157      rot=numpy.dot(refmat(q, -p), refmat(p, -p)) 
158      return rot 
159   
160   
161 -def calc_angle(v1, v2, v3): 
162      """ 
163      Calculate the angle between 3 vectors 
164      representing 3 connected points. 
165   
166      @param v1, v2, v3: the tree points that define the angle 
167      @type v1, v2, v3: L{Vector} 
168   
169      @return: angle 
170      @rtype: float 
171      """ 
172      v1=v1-v2 
173      v3=v3-v2 
174      return v1.angle(v3) 
175   
176   
177 -def calc_dihedral(v1, v2, v3, v4): 
178      """ 
179      Calculate the dihedral angle between 4 vectors 
180      representing 4 connected points. The angle is in 
181      ]-pi, pi]. 
182   
183      @param v1, v2, v3, v4: the four points that define the dihedral angle 
184      @type v1, v2, v3, v4: L{Vector} 
185      """ 
186      ab=v1-v2 
187      cb=v3-v2 
188      db=v4-v3 
189      u=ab**cb 
190      v=db**cb 
191      w=u**v 
192      angle=u.angle(v) 
193      # Determine sign of angle 
194      try: 
195          if cb.angle(w)>0.001: 
196              angle=-angle 
197      except ZeroDivisionError: 
198          # dihedral=pi 
199          pass 
200      return angle 
201   
202   
203 -class Vector(object): 
204      "3D vector" 
205   
206 -    def __init__(self, x, y=None, z=None): 
207          if y is None and z is None: 
208              # Array, list, tuple... 
209              if len(x)!=3: 
210                  raise ValueError("Vector: x is not a " 
211                                   "list/tuple/array of 3 numbers") 
212              self._ar=numpy.array(x, 'd') 
213          else: 
214              # Three numbers 
215              self._ar=numpy.array((x, y, z), 'd') 
216   
217 -    def __repr__(self): 
218          x,y,z=self._ar 
219          return "<Vector %.2f, %.2f, %.2f>" % (x,y,z) 
220   
221 -    def __neg__(self): 
222          "Return Vector(-x, -y, -z)" 
223          a=-self._ar 
224          return Vector(a) 
225   
226 -    def __add__(self, other): 
227          "Return Vector+other Vector or scalar" 
228          if isinstance(other, Vector): 
229              a=self._ar+other._ar 
230          else: 
231              a=self._ar+numpy.array(other) 
232          return Vector(a) 
233   
234 -    def __sub__(self, other): 
235          "Return Vector-other Vector or scalar" 
236          if isinstance(other, Vector): 
237              a=self._ar-other._ar 
238          else: 
239              a=self._ar-numpy.array(other) 
240          return Vector(a) 
241   
242 -    def __mul__(self, other): 
243          "Return Vector.Vector (dot product)" 
244          return sum(self._ar*other._ar) 
245   
246 -    def __div__(self, x): 
247          "Return Vector(coords/a)" 
248          a=self._ar/numpy.array(x) 
249          return Vector(a) 
250   
251 -    def __pow__(self, other): 
252          "Return VectorxVector (cross product) or Vectorxscalar" 
253          if isinstance(other, Vector): 
254              a,b,c=self._ar 
255              d,e,f=other._ar 
256              c1=numpy.linalg.det(numpy.array(((b,c), (e,f)))) 
257              c2=-numpy.linalg.det(numpy.array(((a,c), (d,f)))) 
258              c3=numpy.linalg.det(numpy.array(((a,b), (d,e)))) 
259              return Vector(c1,c2,c3) 
260          else: 
261              a=self._ar*numpy.array(other) 
262              return Vector(a) 
263   
264 -    def __getitem__(self, i): 
265          return self._ar[i] 
266   
267 -    def __setitem__(self, i, value): 
268          self._ar[i]=value 
269   
270 -    def __contains__(self, i): 
271          return (i in self._ar) 
272   
273 -    def norm(self): 
274          "Return vector norm" 
275          return numpy.sqrt(sum(self._ar*self._ar)) 
276   
277 -    def normsq(self): 
278          "Return square of vector norm" 
279          return abs(sum(self._ar*self._ar)) 
280   
281 -    def normalize(self): 
282          "Normalize the Vector" 
283          self._ar=self._ar/self.norm() 
284   
285 -    def normalized(self): 
286          "Return a normalized copy of the Vector" 
287          v=self.copy() 
288          v.normalize() 
289          return v 
290   
291 -    def angle(self, other): 
292          "Return angle between two vectors" 
293          n1=self.norm() 
294          n2=other.norm() 
295          c=(self*other)/(n1*n2) 
296          # Take care of roundoff errors 
297          c=min(c,1) 
298          c=max(-1,c) 
299          return numpy.arccos(c) 
300   
301 -    def get_array(self): 
302          "Return (a copy of) the array of coordinates" 
303          return numpy.array(self._ar) 
304   
305 -    def left_multiply(self, matrix): 
306          "Return Vector=Matrix x Vector" 
307          a=numpy.dot(matrix, self._ar) 
308          return Vector(a) 
309   
310 -    def right_multiply(self, matrix): 
311          "Return Vector=Vector x Matrix" 
312          a=numpy.dot(self._ar, matrix) 
313          return Vector(a) 
314   
315 -    def copy(self): 
316          "Return a deep copy of the Vector" 
317          return Vector(self._ar) 
318   
319  if __name__=="__main__": 
320   
321      from numpy.random import random 
322   
323      v1=Vector(0,0,1) 
324      v2=Vector(0,0,0) 
325      v3=Vector(0,1,0) 
326      v4=Vector(1,1,0) 
327   
328      v4.normalize() 
329   
330      print v4 
331   
332      print calc_angle(v1, v2, v3) 
333      dih=calc_dihedral(v1, v2, v3, v4) 
334      # Test dihedral sign 
335      assert(dih>0) 
336      print "DIHEDRAL ", dih 
337   
338      ref=refmat(v1, v3) 
339      rot=rotmat(v1, v3) 
340   
341      print v3 
342      print v1.left_multiply(ref) 
343      print v1.left_multiply(rot) 
344      print v1.right_multiply(numpy.transpose(rot)) 
345   
346      # - 
347      print v1-v2 
348      print v1-1 
349      print v1+(1,2,3) 
350      # + 
351      print v1+v2 
352      print v1+3 
353      print v1-(1,2,3) 
354      # * 
355      print v1*v2 
356      # / 
357      print v1/2 
358      print v1/(1,2,3) 
359      # ** 
360      print v1**v2 
361      print v1**2 
362      print v1**(1,2,3) 
363      # norm 
364      print v1.norm() 
365      # norm squared 
366      print v1.normsq() 
367      # setitem 
368      v1[2]=10 
369      print v1 
370      # getitem 
371      print v1[2] 
372   
373      print numpy.array(v1) 
374   
375      print "ROT" 
376   
377      angle=random()*numpy.pi 
378      axis=Vector(random(3)-random(3)) 
379      axis.normalize() 
380   
381      m=rotaxis(angle, axis) 
382   
383      cangle, caxis=m2rotaxis(m) 
384   
385      print angle-cangle 
386      print axis-caxis 
387      print 
388