Back to index

python-biopython  1.60
SVDSuperimposer.py
Go to the documentation of this file.
00001 # Copyright (C) 2002, Thomas Hamelryck (thamelry@vub.ac.be)
00002 # This code is part of the Biopython distribution and governed by its
00003 # license.  Please see the LICENSE file that should have been included
00004 # as part of this package.
00005 
00006 from numpy import dot, transpose, sqrt, array
00007 from numpy.linalg import svd, det
00008 
00009 class SVDSuperimposer(object):
00010     """
00011     SVDSuperimposer finds the best rotation and translation to put
00012     two point sets on top of each other (minimizing the RMSD). This is 
00013     eg. useful to superimpose crystal structures.  
00014 
00015     SVD stands for Singular Value Decomposition, which is used to calculate
00016     the superposition.
00017 
00018     Reference:
00019 
00020     Matrix computations, 2nd ed. Golub, G. & Van Loan, CF., The Johns 
00021     Hopkins University Press, Baltimore, 1989
00022     """
00023     def __init__(self):
00024         self._clear()
00025 
00026     # Private methods
00027 
00028     def _clear(self):
00029         self.reference_coords=None
00030         self.coords=None
00031         self.transformed_coords=None
00032         self.rot=None
00033         self.tran=None
00034         self.rms=None
00035         self.init_rms=None
00036 
00037     def _rms(self, coords1, coords2):
00038         "Return rms deviations between coords1 and coords2."
00039         diff=coords1-coords2
00040         l=coords1.shape[0]
00041         return sqrt(sum(sum(diff*diff))/l)
00042 
00043     # Public methods
00044     
00045     def set(self, reference_coords, coords):
00046         """
00047         Set the coordinates to be superimposed.
00048         coords will be put on top of reference_coords.
00049 
00050         o reference_coords: an NxDIM array
00051         o coords: an NxDIM array
00052 
00053         DIM is the dimension of the points, N is the number
00054         of points to be superimposed.
00055         """
00056         # clear everything from previous runs
00057         self._clear()
00058         # store cordinates
00059         self.reference_coords=reference_coords
00060         self.coords=coords
00061         n=reference_coords.shape
00062         m=coords.shape
00063         if n!=m or not(n[1]==m[1]==3):
00064             raise Exception("Coordinate number/dimension mismatch.")
00065         self.n=n[0]
00066 
00067     def run(self):
00068         "Superimpose the coordinate sets."
00069         if self.coords is None or self.reference_coords is None:
00070             raise Exception("No coordinates set.")
00071         coords=self.coords
00072         reference_coords=self.reference_coords
00073         # center on centroid
00074         av1=sum(coords)/self.n  
00075         av2=sum(reference_coords)/self.n    
00076         coords=coords-av1
00077         reference_coords=reference_coords-av2
00078         # correlation matrix
00079         a=dot(transpose(coords), reference_coords)
00080         u, d, vt=svd(a)
00081         self.rot=transpose(dot(transpose(vt), transpose(u)))
00082         # check if we have found a reflection
00083         if det(self.rot)<0:
00084             vt[2]=-vt[2]
00085             self.rot=transpose(dot(transpose(vt), transpose(u)))
00086         self.tran=av2-dot(av1, self.rot)
00087 
00088     def get_transformed(self):
00089         "Get the transformed coordinate set."
00090         if self.coords is None or self.reference_coords is None:
00091             raise Exception("No coordinates set.")
00092         if self.rot is None:
00093             raise Exception("Nothing superimposed yet.")
00094         if self.transformed_coords is None:
00095             self.transformed_coords=dot(self.coords, self.rot)+self.tran
00096         return self.transformed_coords
00097 
00098     def get_rotran(self):
00099         "Right multiplying rotation matrix and translation."
00100         if self.rot is None:
00101             raise Exception("Nothing superimposed yet.")
00102         return self.rot, self.tran
00103 
00104     def get_init_rms(self):
00105         "Root mean square deviation of untransformed coordinates."
00106         if self.coords is None:
00107             raise Exception("No coordinates set yet.")
00108         if self.init_rms is None:
00109             self.init_rms=self._rms(self.coords, self.reference_coords)
00110         return self.init_rms
00111 
00112     def get_rms(self):
00113         "Root mean square deviation of superimposed coordinates."
00114         if self.rms is None:
00115             transformed_coords=self.get_transformed()
00116             self.rms=self._rms(transformed_coords, self.reference_coords)
00117         return self.rms
00118 
00119 
00120 if __name__=="__main__":
00121 
00122     # start with two coordinate sets (Nx3 arrays - float)
00123 
00124     x=array([[51.65, -1.90, 50.07],
00125          [50.40, -1.23, 50.65],
00126          [50.68, -0.04, 51.54],
00127          [50.22, -0.02, 52.85]], 'f')
00128 
00129     y=array([[51.30, -2.99, 46.54],
00130          [51.09, -1.88, 47.58],
00131          [52.36, -1.20, 48.03],
00132          [52.71, -1.18, 49.38]], 'f')
00133 
00134     # start!
00135     sup=SVDSuperimposer()
00136 
00137     # set the coords
00138     # y will be rotated and translated on x
00139     sup.set(x, y)
00140 
00141     # do the lsq fit
00142     sup.run()
00143 
00144     # get the rmsd
00145     rms=sup.get_rms()
00146 
00147     # get rotation (right multiplying!) and the translation
00148     rot, tran=sup.get_rotran()
00149 
00150     # rotate y on x
00151     y_on_x1=dot(y, rot)+tran
00152 
00153     # same thing
00154     y_on_x2=sup.get_transformed()
00155 
00156     print y_on_x1
00157     print
00158     print y_on_x2
00159     print
00160     print "%.2f" % rms