Package Bio :: Package PDB :: Module Vector'
[hide private]
[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 
  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  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 92 theta rad around vector. 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
245 - def __add__(self, other):
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