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 # 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 """Vector to axis method. 56 57 Return the vector between a point and 58 the closest point on a line (ie. the perpendicular 59 projection of the point on the line). 60 61 @type line: L{Vector} 62 @param line: vector defining a line 63 64 @type point: L{Vector} 65 @param point: vector defining the point 66 """ 67 line = line.normalized() 68 np = point.norm() 69 angle = line.angle(point) 70 return point - line ** (np * numpy.cos(angle))
71 72
73 -def rotaxis2m(theta, vector):
74 """Calculate left multiplying rotation matrix. 75 76 Calculate a left multiplying rotation matrix that rotates 77 theta rad around vector. 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 115 rotaxis = rotaxis2m 116 117
118 -def refmat(p, q):
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 """Calculate angle method. 167 168 Calculate the angle between 3 vectors 169 representing 3 connected points. 170 171 @param v1, v2, v3: the tree points that define the angle 172 @type v1, v2, v3: L{Vector} 173 174 @return: angle 175 @rtype: float 176 """ 177 v1 = v1 - v2 178 v3 = v3 - v2 179 return v1.angle(v3)
180 181
182 -def calc_dihedral(v1, v2, v3, v4):
183 """Calculate dihedral angle method. 184 185 Calculate the dihedral angle between 4 vectors 186 representing 4 connected points. The angle is in 187 ]-pi, pi]. 188 189 @param v1, v2, v3, v4: the four points that define the dihedral angle 190 @type v1, v2, v3, v4: L{Vector} 191 """ 192 ab = v1 - v2 193 cb = v3 - v2 194 db = v4 - v3 195 u = ab ** cb 196 v = db ** cb 197 w = u ** v 198 angle = u.angle(v) 199 # Determine sign of angle 200 try: 201 if cb.angle(w) > 0.001: 202 angle = -angle 203 except ZeroDivisionError: 204 # dihedral=pi 205 pass 206 return angle
207 208
209 -class Vector(object):
210 """3D vector.""" 211
212 - def __init__(self, x, y=None, z=None):
213 if y is None and z is None: 214 # Array, list, tuple... 215 if len(x) != 3: 216 raise ValueError("Vector: x is not a " 217 "list/tuple/array of 3 numbers") 218 self._ar = numpy.array(x, 'd') 219 else: 220 # Three numbers 221 self._ar = numpy.array((x, y, z), 'd')
222
223 - def __repr__(self):
224 x, y, z = self._ar 225 return "<Vector %.2f, %.2f, %.2f>" % (x, y, z)
226
227 - def __neg__(self):
228 """Return Vector(-x, -y, -z).""" 229 a = -self._ar 230 return Vector(a)
231
232 - def __add__(self, other):
233 """Return Vector+other Vector or scalar.""" 234 if isinstance(other, Vector): 235 a = self._ar + other._ar 236 else: 237 a = self._ar + numpy.array(other) 238 return Vector(a)
239
240 - def __sub__(self, other):
241 """Return Vector-other Vector or scalar.""" 242 if isinstance(other, Vector): 243 a = self._ar - other._ar 244 else: 245 a = self._ar - numpy.array(other) 246 return Vector(a)
247
248 - def __mul__(self, other):
249 """Return Vector.Vector (dot product).""" 250 return sum(self._ar * other._ar)
251
252 - def __div__(self, x):
253 """Return Vector(coords/a).""" 254 a = self._ar / numpy.array(x) 255 return Vector(a)
256
257 - def __pow__(self, other):
258 """Return VectorxVector (cross product) or Vectorxscalar.""" 259 if isinstance(other, Vector): 260 a, b, c = self._ar 261 d, e, f = other._ar 262 c1 = numpy.linalg.det(numpy.array(((b, c), (e, f)))) 263 c2 = -numpy.linalg.det(numpy.array(((a, c), (d, f)))) 264 c3 = numpy.linalg.det(numpy.array(((a, b), (d, e)))) 265 return Vector(c1, c2, c3) 266 else: 267 a = self._ar * numpy.array(other) 268 return Vector(a)
269
270 - def __getitem__(self, i):
271 return self._ar[i]
272
273 - def __setitem__(self, i, value):
274 self._ar[i] = value
275
276 - def __contains__(self, i):
277 return (i in self._ar)
278
279 - def norm(self):
280 """Return vector norm.""" 281 return numpy.sqrt(sum(self._ar * self._ar))
282
283 - def normsq(self):
284 """Return square of vector norm.""" 285 return abs(sum(self._ar * self._ar))
286
287 - def normalize(self):
288 """Normalize the Vector.""" 289 self._ar = self._ar / self.norm()
290
291 - def normalized(self):
292 """Return a normalized copy of the Vector.""" 293 v = self.copy() 294 v.normalize() 295 return v
296
297 - def angle(self, other):
298 """Return angle between two vectors.""" 299 n1 = self.norm() 300 n2 = other.norm() 301 c = (self * other) / (n1 * n2) 302 # Take care of roundoff errors 303 c = min(c, 1) 304 c = max(-1, c) 305 return numpy.arccos(c)
306
307 - def get_array(self):
308 """Return (a copy of) the array of coordinates.""" 309 return numpy.array(self._ar)
310
311 - def left_multiply(self, matrix):
312 """Return Vector=Matrix x Vector.""" 313 a = numpy.dot(matrix, self._ar) 314 return Vector(a)
315
316 - def right_multiply(self, matrix):
317 """Return Vector=Vector x Matrix.""" 318 a = numpy.dot(self._ar, matrix) 319 return Vector(a)
320
321 - def copy(self):
322 """Return a deep copy of the Vector.""" 323 return Vector(self._ar)
324