[frames] | no frames]

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

```

 Generated by Epydoc 3.0.1 on Thu Aug 25 13:17:00 2016 http://epydoc.sourceforge.net