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 """ 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 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 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
230 - def __add__(self, other):
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