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

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