Package Bio :: Package SVDSuperimposer
[hide private]
[frames] | no frames]

Source Code for Package Bio.SVDSuperimposer

  1  # Copyright (C) 2002, Thomas Hamelryck (thamelry@vub.ac.be) 
  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  SVDSuperimposer finds the best rotation and translation to put 
  7  two point sets on top of each other (minimizing the RMSD). This is 
  8  eg. useful to superimpose crystal structures. SVD stands for singular 
  9  value decomposition, which is used in the algorithm. 
 10  """ 
 11   
 12  from __future__ import print_function 
 13   
 14  from numpy import dot, transpose, sqrt, array 
 15  from numpy.linalg import svd, det 
 16   
 17   
18 -class SVDSuperimposer(object):
19 """ 20 SVDSuperimposer finds the best rotation and translation to put 21 two point sets on top of each other (minimizing the RMSD). This is 22 eg. useful to superimpose crystal structures. 23 24 SVD stands for Singular Value Decomposition, which is used to calculate 25 the superposition. 26 27 Reference: 28 29 Matrix computations, 2nd ed. Golub, G. & Van Loan, CF., The Johns 30 Hopkins University Press, Baltimore, 1989 31 """
32 - def __init__(self):
33 self._clear()
34 35 # Private methods 36
37 - def _clear(self):
38 self.reference_coords = None 39 self.coords = None 40 self.transformed_coords = None 41 self.rot = None 42 self.tran = None 43 self.rms = None 44 self.init_rms = None
45
46 - def _rms(self, coords1, coords2):
47 "Return rms deviations between coords1 and coords2." 48 diff = coords1 - coords2 49 l = coords1.shape[0] 50 return sqrt(sum(sum(diff*diff))/l)
51 52 # Public methods 53
54 - def set(self, reference_coords, coords):
55 """ 56 Set the coordinates to be superimposed. 57 coords will be put on top of reference_coords. 58 59 o reference_coords: an NxDIM array 60 o coords: an NxDIM array 61 62 DIM is the dimension of the points, N is the number 63 of points to be superimposed. 64 """ 65 # clear everything from previous runs 66 self._clear() 67 # store cordinates 68 self.reference_coords = reference_coords 69 self.coords = coords 70 n = reference_coords.shape 71 m = coords.shape 72 if n != m or not(n[1] == m[1] == 3): 73 raise Exception("Coordinate number/dimension mismatch.") 74 self.n = n[0]
75
76 - def run(self):
77 "Superimpose the coordinate sets." 78 if self.coords is None or self.reference_coords is None: 79 raise Exception("No coordinates set.") 80 coords = self.coords 81 reference_coords = self.reference_coords 82 # center on centroid 83 av1 = sum(coords) / self.n 84 av2 = sum(reference_coords) / self.n 85 coords = coords - av1 86 reference_coords = reference_coords - av2 87 # correlation matrix 88 a = dot(transpose(coords), reference_coords) 89 u, d, vt = svd(a) 90 self.rot = transpose(dot(transpose(vt), transpose(u))) 91 # check if we have found a reflection 92 if det(self.rot) < 0: 93 vt[2] = -vt[2] 94 self.rot = transpose(dot(transpose(vt), transpose(u))) 95 self.tran = av2 - dot(av1, self.rot)
96
97 - def get_transformed(self):
98 "Get the transformed coordinate set." 99 if self.coords is None or self.reference_coords is None: 100 raise Exception("No coordinates set.") 101 if self.rot is None: 102 raise Exception("Nothing superimposed yet.") 103 if self.transformed_coords is None: 104 self.transformed_coords = dot(self.coords, self.rot) + self.tran 105 return self.transformed_coords
106
107 - def get_rotran(self):
108 "Right multiplying rotation matrix and translation." 109 if self.rot is None: 110 raise Exception("Nothing superimposed yet.") 111 return self.rot, self.tran
112
113 - def get_init_rms(self):
114 "Root mean square deviation of untransformed coordinates." 115 if self.coords is None: 116 raise Exception("No coordinates set yet.") 117 if self.init_rms is None: 118 self.init_rms = self._rms(self.coords, self.reference_coords) 119 return self.init_rms
120
121 - def get_rms(self):
122 "Root mean square deviation of superimposed coordinates." 123 if self.rms is None: 124 transformed_coords = self.get_transformed() 125 self.rms = self._rms(transformed_coords, self.reference_coords) 126 return self.rms
127 128 129 if __name__ == "__main__": 130 131 # start with two coordinate sets (Nx3 arrays - float) 132 133 x = array([[51.65, -1.90, 50.07], 134 [50.40, -1.23, 50.65], 135 [50.68, -0.04, 51.54], 136 [50.22, -0.02, 52.85]], 'f') 137 138 y = array([[51.30, -2.99, 46.54], 139 [51.09, -1.88, 47.58], 140 [52.36, -1.20, 48.03], 141 [52.71, -1.18, 49.38]], 'f') 142 143 # start! 144 sup = SVDSuperimposer() 145 146 # set the coords 147 # y will be rotated and translated on x 148 sup.set(x, y) 149 150 # do the lsq fit 151 sup.run() 152 153 # get the rmsd 154 rms = sup.get_rms() 155 156 # get rotation (right multiplying!) and the translation 157 rot, tran = sup.get_rotran() 158 159 # rotate y on x 160 y_on_x1 = dot(y, rot) + tran 161 162 # same thing 163 y_on_x2 = sup.get_transformed() 164 165 print(y_on_x1) 166 print("") 167 print(y_on_x2) 168 print("") 169 print("%.2f" % rms) 170