[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
4  # as part of this package.
5
6  from numpy import dot, transpose, sqrt, array
7  from numpy.linalg import svd, det
8
9
10 -class SVDSuperimposer(object):
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 = dot(transpose(coords), reference_coords)
81          u, d, vt = svd(a)
82          self.rot = transpose(dot(transpose(vt), transpose(u)))
83          # check if we have found a reflection
84          if det(self.rot) < 0:
85              vt[2] = -vt[2]
86              self.rot = transpose(dot(transpose(vt), transpose(u)))
87          self.tran = av2 - dot(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 = dot(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
124
125      x = array([[51.65, -1.90, 50.07],
126           [50.40, -1.23, 50.65],
127           [50.68, -0.04, 51.54],
128           [50.22, -0.02, 52.85]], 'f')
129
130      y = array([[51.30, -2.99, 46.54],
131           [51.09, -1.88, 47.58],
132           [52.36, -1.20, 48.03],
133           [52.71, -1.18, 49.38]], 'f')
134
135      # start!
136      sup = SVDSuperimposer()
137
138      # set the coords
139      # y will be rotated and translated on x
140      sup.set(x, y)
141
142      # do the lsq fit
143      sup.run()
144
145      # get the rmsd
146      rms = sup.get_rms()
147
148      # get rotation (right multiplying!) and the translation
149      rot, tran = sup.get_rotran()
150
151      # rotate y on x
152      y_on_x1 = dot(y, rot) + tran
153
154      # same thing
155      y_on_x2 = sup.get_transformed()
156
157      print y_on_x1
158      print
159      print y_on_x2
160      print
161      print "%.2f" % rms
162
<!--
expandto(location.href);
// -->

```

 Generated by Epydoc 3.0.1 on Tue Feb 5 18:03:10 2013 http://epydoc.sourceforge.net