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

```

 Generated by Epydoc 3.0.1 on Mon Jul 10 15:16:51 2017 http://epydoc.sourceforge.net