Source Code for Module Biskit.rmsFit

  1  ## 
  2  ## Biskit, a toolkit for the manipulation of macromolecular structures 
  3  ## Copyright (C) 2002-2005 Michael Habeck 
  4  ## 
  5  ## This program is free software; you can redistribute it and/or 
  6  ## modify it under the terms of the GNU General Public License as 
  7  ## published by the Free Software Foundation; either version 2 of the 
  8  ## License, or any later version. 
  9  ## 
 10  ## This program is distributed in the hope that it will be useful, 
 11  ## but WITHOUT ANY WARRANTY; without even the implied warranty of 
 12  ## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU 
 13  ## General Public License for more details. 
 14  ## 
 15  ## You find a copy of the GNU General Public License in the file 
 16  ## license.txt along with this program; if not, write to the Free 
 17  ## Software Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA. 
 18   
 19  ## Contributions: Michael Habeck, Raik Gruenberg, Johan Leckner 
 20  ## last $Author: leckner $ 
 21  ## last $Date: 2006/09/12 06:34:36 $ 
 22  ## $Revision: 2.5 $ 
 23  """ 
 24  superimpose 2 structures iteratively 
 25  """ 
 26   
 27  import Biskit.mathUtils as MU 
 28  import Numeric as N 
 29  from LinearAlgebra import singular_value_decomposition as svd 
 30   
 31  ## def average(x): 
 32  ##     return N.sum(x) / len(x) 
 33   
 34  ## def variance(x): 
 35  ##     return N.average(N.power(x - N.average(x), 2)) 
 36   
 37  ## def standardDev(x): 
 38  ##     return N.sqrt(variance(x)) 
 39   
 40   
 41 -def findTransformation(x, y): 
 42      """ 
 43      Match two arrays by rotation and translation. Returns the 
 44      rotation matrix and the translation vector. 
 45   
 46      @param x: first set of coordinates 
 47      @type  x: array('f') 
 48      @param y: second set of coordinates 
 49      @type  y: array('f') 
 50       
 51      @return: rotation matrix (3x3) and translation vector (1x3) 
 52      @rtype:  array, array 
 53      """ 
 54      ## center configurations 
 55      x_av = N.average(x) 
 56      y_av = N.average(y) 
 57   
 58      x = x - x_av 
 59      y = y - y_av 
 60   
 61      ## svd of correlation matrix 
 62      v, l, u = svd(N.dot(N.transpose(x), y)) 
 63   
 64      ## build rotation matrix and translation vector 
 65      r = N.dot(v, u) 
 66      t = x_av - N.dot(r, y_av) 
 67   
 68      return r, t 
 69   
 70   
 71 -def match(x, y, n_iterations=1, z=2, eps_rmsd=0.5, eps_stdv=0.05): 
 72      """ 
 73      Matches two arrays onto each other, while iteratively removing outliers. 
 74      Superimposed array y would be C{ N.dot(y, N.transpose(r)) + t }. 
 75       
 76      @param n_iterations: number of calculations:: 
 77                             1 .. no iteration  
 78                             0 .. until convergence 
 79      @type  n_iterations: 1|0 
 80      @param z: number of standard deviations for outlier definition (default: 2) 
 81      @type  z: float 
 82      @param eps_rmsd: tolerance in rmsd (default: 0.5) 
 83      @type  eps_rmsd: float 
 84      @param eps_stdv: tolerance in standard deviations (default: 0.05) 
 85      @type  eps_stdv: float 
 86       
 87      @return: (r,t), [ [percent_considered, rmsd_for_it, outliers] ] 
 88      @rtype: (array, array), [float, float, int] 
 89      """ 
 90      iter_trace = [] 
 91   
 92      rmsd_old = 0 
 93      stdv_old = 0 
 94   
 95      n = 0 
 96      converged = 0 
 97   
 98      mask = N.ones(len(y)) 
 99   
100      while not converged: 
101   
102          ## find transformation for best match 
103          r, t = findTransformation(N.compress(mask, x, 0), 
104                                    N.compress(mask, y, 0)) 
105   
106          ## transform coordinates 
107          xt = N.dot(y, N.transpose(r)) + t 
108   
109          ## calculate row distances 
110          d = N.sqrt(N.sum(N.power(x - xt, 2), 1)) * mask 
111   
112          ## calculate rmsd and stdv 
113          rmsd = N.sqrt(N.average(N.compress(mask, d)**2)) 
114          stdv = MU.SD(N.compress(mask, d)) 
115   
116          ## check conditions for convergence 
117          d_rmsd = abs(rmsd - rmsd_old) 
118          d_stdv = abs(1 - stdv_old / stdv) 
119   
120          if d_rmsd < eps_rmsd and d_stdv < eps_stdv: 
121              converged = 1 
122          else: 
123              rmsd_old = rmsd 
124              stdv_old = stdv 
125   
126          ## store result 
127          perc = round(float(N.sum(mask)) / float(len(mask)), 2) 
128   
129          ## throw out non-matching rows 
130          mask = N.logical_and(mask, N.less(d, rmsd + z * stdv)) 
131          outliers = N.nonzero( N.logical_not( mask ) ) 
132          iter_trace.append([perc, round(rmsd, 3), outliers]) 
133   
134          n += 1 
135   
136          if n_iterations and n >= n_iterations: 
137              break 
138   
139      return (r, t), iter_trace 
140   
141   
142 -def rowDistances( x, y ): 
143      """ 
144      Calculate the distances between the items of two arrays (of same shape) 
145      after least-squares superpositioning. 
146       
147      @param x: first set of coordinates 
148      @type  x: array('f') 
149      @param y: second set of coordinates 
150      @type  y: array('f')   
151       
152      @return: array( len(x), 'f' ), distance between x[i] and y[i] for all i 
153      @rtype: array 
154      """ 
155      ## find transformation for best match 
156      r, t = findTransformation(x, y) 
157   
158      ## transform coordinates 
159      z = N.dot(y, N.transpose(r)) + t 
160   
161      ## calculate row distances 
162      return N.sqrt(N.sum(N.power(x - z, 2), 1))  
163   
164   
165   
166  ############# 
167  ##  TESTING         
168  ############# 
169           
170 -class Test: 
171      """ 
172      Test class 
173      """ 
174   
175       
176 -    def run( self, local=0 ): 
177          """ 
178          run function test 
179           
180          @param local: transfer local variables to global and perform 
181                        other tasks only when run locally 
182          @type  local: 1|0 
183           
184          @return: rotation matrix 
185          @rtype: array 
186          """ 
187          import Biskit.tools as T 
188   
189          traj = T.Load( T.testRoot() + '/lig_pcr_00/traj.dat' ) 
190   
191          rt, rmsdLst = match( traj.ref.xyz, traj[-1].xyz) 
192   
193          if local: 
194              print 'RMSD: %.2f'%rmsdLst[0][1] 
195              globals().update( locals() ) 
196           
197          # return rotation matrix 
198          return rt[0] 
199   
200   
201 -    def expected_result( self ): 
202          """ 
203          Precalculated result to check for consistent performance. 
204   
205          @return: rotation matrix 
206          @rtype:  array 
207          """ 
208          return N.array( [[ 0.9999011,   0.01311352,  0.00508244,], 
209                           [-0.01310219,  0.99991162, -0.00225578,], 
210                           [-0.00511157,  0.00218896,  0.99998454 ]] ) 
211           
212   
213  if __name__ == '__main__': 
214   
215      test = Test() 
216   
217      assert abs( N.sum( N.ravel( test.run( local=1 )- test.expected_result() ) ) ) < 1e-6 
218