IMP  2.0.1
The Integrative Modeling Platform
distance.h
Go to the documentation of this file.
1 /**
2  * \file IMP/atom/distance.h
3  * \brief distance metrics
4  *
5  * Copyright 2007-2013 IMP Inventors. All rights reserved.
6  *
7  */
8 
9 #ifndef IMPATOM_DISTANCE_H
10 #define IMPATOM_DISTANCE_H
11 
12 #include <IMP/atom/atom_config.h>
13 #include <IMP/core/XYZ.h>
14 #include "Hierarchy.h"
15 #include "Selection.h"
16 #include "IMP/base_types.h"
17 
18 IMPATOM_BEGIN_NAMESPACE
19 
20 //! Calculate the root mean square deviation between two sets of 3D points.
21 /**
22  \note the function assumes correspondence between the two sets of
23  points and does not perform rigid alignment.
24 
25  \genericgeometry
26  */
27 template <class Vector3DsOrXYZs0, class Vector3DsOrXYZs1>
28 inline double get_rmsd(const Vector3DsOrXYZs0& m1, const Vector3DsOrXYZs1& m2,
29  const IMP::algebra::Transformation3D &tr_for_second
30  = IMP::algebra::get_identity_transformation_3d()) {
31  IMP_USAGE_CHECK(std::distance(m1.begin(), m1.end())
32  ==std::distance(m2.begin(), m2.end()),
33  "The input sets of XYZ points "
34  <<"should be of the same size");
35  double rmsd=0.0;
36  typename Vector3DsOrXYZs0::const_iterator it0= m1.begin();
37  typename Vector3DsOrXYZs1::const_iterator it1= m2.begin();
38  for(; it0!= m1.end(); ++it0, ++it1) {
39  algebra::Vector3D tred
40  =tr_for_second.get_transformed(get_vector_d_geometry(*it1));
41  rmsd += algebra::get_squared_distance(get_vector_d_geometry(*it0),
42  tred);
43  }
44  return std::sqrt(rmsd / m1.size());
45 }
46 
47 /** RMSD on a pair of Selections.*/
48 inline double get_rmsd(const Selection &s0,
49  const Selection &s1,
50  const IMP::algebra::Transformation3D &tr_for_second
51  = IMP::algebra::get_identity_transformation_3d()) {
52  return get_rmsd(s0.get_selected_particles(),
53  s1.get_selected_particles(), tr_for_second);
54 }
55 
56 //! Calculate the root mean square deviation between two sets of 3D points.
57 /**
58  \note the function assumes correspondence between the two sets of
59  points and does not need rigid alignment. Note that it is different from
60  get_drms().
61 
62  \genericgeometry
63  */
64 template <class Vector3DsOrXYZs0, class Vector3DsOrXYZs1>
65 inline double get_drmsd(const Vector3DsOrXYZs0& m0, const Vector3DsOrXYZs1& m1)
66  {
67  IMP_USAGE_CHECK(m0.size()==m1.size(),
68  "The input sets of XYZ points "
69  <<"should be of the same size");
70  double drmsd=0.0;
71 
72  int npairs = 0;
73  for(unsigned i=0;i<m0.size()-1;++i){
74  algebra::Vector3D v0i=get_vector_d_geometry(m0[i]);
75  algebra::Vector3D v1i=get_vector_d_geometry(m1[i]);
76 
77  for(unsigned j=i+1;j<m0.size();++j){
78  algebra::Vector3D v0j=get_vector_d_geometry(m0[j]);
79  algebra::Vector3D v1j=get_vector_d_geometry(m1[j]);
80 
81  double dist0=algebra::get_distance(v0i,v0j);
82  double dist1=algebra::get_distance(v1i,v1j);
83  drmsd+=(dist0-dist1)*(dist0-dist1);
84  npairs++;
85  }
86  }
87  return std::sqrt(drmsd/npairs);
88 }
89 
90 
91 //! Computes the native overlap between two sets of 3D points
92 /**
93  \param[in] m1 first set
94  \param[in] m2 second set
95  \param[in] threshold threshold distance (angstroms) for the calculation
96  \note the function assumes correspondence between two sets of points and does
97  not perform rigid alignment.
98  \return the native overlap, as a percentage (from 0 to 100)
99  \genericgeometry
100 **/
101 template <class Vector3DsOrXYZs0, class Vector3DsOrXYZs1>
102 inline double get_native_overlap(const Vector3DsOrXYZs0& m1,
103  const Vector3DsOrXYZs1& m2,double threshold) {
104  IMP_USAGE_CHECK(m1.size()==m2.size(),
105  "native_verlap: The input sets of XYZ points "
106  <<"should be of the same size");
107  unsigned int distances=0;
108  for(unsigned int i=0;i<m1.size();i++) {
109  double d = algebra::get_distance(get_vector_d_geometry(m1[i]),
110  get_vector_d_geometry(m2[i]));
111  if(d<=threshold) distances++;
112  }
113  return 100.0*distances/m1.size();
114 }
115 
116 
117 /*! Distance-RMS between two sets of points, defined as:
118  sqrt( sum[ (d1ij**2 - d2ij**2)**2]/(4 * sum[d1ij**2]) )
119  (Levitt, 1992)
120  d1ij - distance between points i and j in set m1 (the "reference" set)
121  d2ij - distance between points i and j in set m2
122  \param[in] m1 set of points
123  \param[in] m2 set of points
124 */
125 template <class Vector3DsOrXYZs0, class Vector3DsOrXYZs1>
126 inline double get_drms(const Vector3DsOrXYZs0& m1,
127  const Vector3DsOrXYZs1& m2) {
128  IMP_USAGE_CHECK(m1.size()==m2.size(),
129  "native_verlap: The input sets of XYZ points "
130  <<"should be of the same size");
131 
132  unsigned int n = m1.size();
133  double drms = 0.0;
134  double sum_d1ij = 0.0;
135  for (unsigned int i=0; i < n; ++i) {
136  for (unsigned int j= i+1; j< n; ++j) {
137  double sqd1 = algebra::get_squared_distance(get_vector_d_geometry(m1[i]),
138  get_vector_d_geometry(m1[j]));
139  double sqd2 = algebra::get_squared_distance(get_vector_d_geometry(m2[i]),
140  get_vector_d_geometry(m2[j]));
141  drms += (sqd1 - sqd2)*(sqd1 - sqd2);
142  sum_d1ij += sqd1;
143  }
144  }
145  drms /= (4 * sum_d1ij);
146  return sqrt(drms);
147 }
148 
149 
150 /*! DRMS between to sets of rigid bodies. Points ij belonging to the same
151 rigid body are not evaluated, as they are the same in both sets
152  \param[in] m1 set of points
153  \param[in] m2 set of points
154  \param[in] ranges of points considered to be the same rigid body. Eg,
155  (0-29,30-49) means two rigid bodies, first with the 30 first atoms,
156  second with the last 20
157 */
158 template <class Vector3DsOrXYZs0, class Vector3DsOrXYZs1>
159 inline double get_rigid_bodies_drms(const Vector3DsOrXYZs0 &m1,
160  const Vector3DsOrXYZs1 &m2,
161  const IMP::IntRanges &ranges) {
162  IMP_USAGE_CHECK(m1.size() == m2.size(),
163  "get_rigid_bodies_drms: The input sets of XYZ points "
164  <<"should be of the same size");
165  int n = m1.size();
166  int rn = ranges.size();
167  double drms = 0.0;
168  double sum_d1ij = 0.0;
169  for ( int i=0; i < n; ++i) {
170  int range1 = -1;
171  for ( int k=0; k < rn; ++k) {
172  if(i >= ranges[k].first && i <= ranges[k].second) {
173  range1 = k;
174  break;
175  }
176  }
177  IMP_USAGE_CHECK(range1 >=0, "Point " << i << " of m1 does not belong to "
178  "any range");
179  for ( int j= i+1; j < n; ++j) {
180  int range2 = -1;
181  for ( int k=0; k < rn; ++k) {
182  if(j >= ranges[k].first && j <= ranges[k].second) {
183  range2 = k;
184  break;
185  }
186  }
187  IMP_USAGE_CHECK(range2 >=0, "Point " << j << " of m2 does not belong to "
188  "any range");
189 
190  double sqd1 = algebra::get_squared_distance(get_vector_d_geometry(m1[i]),
191  get_vector_d_geometry(m1[j]));
192  sum_d1ij += sqd1;
193  if(range1 != range2) {
194  // points i and j in different ranges compare distances
195  double sqd2 = algebra::get_squared_distance(
196  get_vector_d_geometry(m2[i]),
197  get_vector_d_geometry(m2[j]));
198  drms += (sqd1 - sqd2) * (sqd1 - sqd2);
199  }
200  }
201  }
202  drms /= (4 * sum_d1ij);
203  return sqrt(drms);
204 }
205 
206 
207 
208 
209 //! Measure the difference between two placements of the same set of points
210 /**
211  \param[in] source The reference placement represented by XYZ coordinates
212  \param[in] target The modeled placement represented by XYZ coordinates
213  \note The measure quantifies the difference between placements
214  of the same structure. A rigid transformation that brings mdl1 to
215  ref1 is reported.
216  \return (d,a), A transformation from mdl to ref represented by
217  a a distance (d) and an angle (a).
218  d is the distance bewteen the centroids of the two
219  placements and a is the axis angle of the rotation matrix between
220  the two placements
221  \note see Lasker,Topf et al JMB, 2009 for details
222  */
223 IMPATOMEXPORT FloatPair
224 get_placement_score(const core::XYZs& source,
225  const core::XYZs& target);
226 
227 //! Measure the difference between two placements of the same set of points
228 /**
229 \param[in] ref1 The reference placement of the first component represented
230  by XYZ coordinates
231 \param[in] ref2 The reference placement of the second component represented
232  by XYZ coordinates
233 \param[in] mdl1 The modeled placement of the first component
234  represented by XYZ coordinates
235 \param[in] mdl2 The modeled placement of the second component
236  represented by XYZ coordinates
237 \return the function returns (distance score,angle score)
238 \note The measure quantifies the difference between the relative placements
239  of two components compared to a reference relative placement.
240  First, the two compared structures are brought into
241  the same frame of reference by superposing the first pair of
242  equivalent domains (ref1 and mdl1). Next, the distance and angle
243  scores are calculated for the second component using placement_score.
244 \note see Topf, Lasker et al Structure, 2008 for details
245  */
246 IMPATOMEXPORT FloatPair get_component_placement_score(
247  const core::XYZs& ref1 ,const core::XYZs& ref2,
248  const core::XYZs& mdl1 ,const core::XYZs& mdl2);
249 
250 //! Measure the RMSD between two placements of the same set of points
251 /**
252 \param[in] ref1 The reference placement of the first component represented
253  by XYZ coordinates
254 \param[in] ref2 The reference placement of the second component represented
255  by XYZ coordinates
256 \param[in] mdl1 The modeled placement of the first component
257  represented by XYZ coordinates
258 \param[in] mdl2 The modeled placement of the second component
259  represented by XYZ coordinates
260 \note The measure quantifies the RMSD between the relative placements
261  of two components compared to a reference relative placement.
262  First, the two compared structures are brought into
263  the same frame of reference by superposing the first pair of
264  equivalent domains (ref1 and mdl1). Next, the RMSD is calculated
265  for the second component
266 \note see Lasker et al JMB, 2009 for details
267  */
268 IMPATOMEXPORT double get_pairwise_rmsd_score(
269  const core::XYZs& ref1 ,const core::XYZs& ref2,
270  const core::XYZs& mdl1 ,const core::XYZs& mdl2);
271 
272 
273 /** Compute the radius of gyration of a set of particles with (optional)
274  radii and mass and (non-optional) coordinates. Either all particles
275  must have mass or none of them.
276 */
277 IMPATOMEXPORT double get_radius_of_gyration(const ParticlesTemp &ps);
278 
279 /** \brief Used to calculate rmsd
280  between multiple transformation that operate on the same particles
281 
282  The constructor takes a list of particles defining the thing which
283  is having its rmsd calculated.
284 */
285 class IMPATOMEXPORT RMSDCalculator {
286 public:
287  RMSDCalculator() {}
288  //! Constructor
289  /**
290  \param[in] ps the particles on which the transformation operate.
291  RMSD will be calculate on these particles.
292  */
293  RMSDCalculator(const ParticlesTemp &ps);
294 
295  //! Get rmsd between two transformations
296  /** The number returned is equilavent to the pseudo code
297 \code
298  get_rmsd(t1.get_transformed(ps),
299  t2.get_transformed(ps));
300 \endcode
301  where ps is the list of particles passed to the constructor.
302  */
303  double get_rmsd(const algebra::Transformation3D& t1,
304  const algebra::Transformation3D& t2) {
305  return sqrt(get_squared_rmsd(t1, t2));}
306 
307  //! Get the squared rmsd between two transformations
308  double get_squared_rmsd(const algebra::Transformation3D& t1,
309  const algebra::Transformation3D& t2);
310  IMP_SHOWABLE_INLINE(RMSDCalculator, out << centroid_);
311 private:
312  algebra::Vector3D centroid_;
313  double d_[3][3];//partial calculation
314  //for example dist[0][1] is the dot product of a two vectors of lenght N
315  //one of all X coordiantes and the second of all Y coordiantes.
316  //N is the number of particles
317 };
318 
319 IMP_VALUES(RMSDCalculator, RMSDCalculators);
320 
321 IMPATOM_END_NAMESPACE
322 
323 #endif /* IMPATOM_DISTANCE_H */