IMP  2.2.1
The Integrative Modeling Platform
atom/distance.h
Go to the documentation of this file.
1 /**
2  * \file IMP/atom/distance.h
3  * \brief distance metrics
4  *
5  * Copyright 2007-2014 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/algebra/distance.h>
17 #include "IMP/base_types.h"
18 
19 IMPATOM_BEGIN_NAMESPACE
20 /**
21  Get the rmsd between two lists of particles
22  */
23 IMPATOMEXPORT double get_rmsd_transforming_first(
24  const IMP::algebra::Transformation3D& tr, const core::XYZs& s0,
25  const core::XYZs& s1);
26 
27 /**
28  Get the rmsd between two lists of particles
29  */
30 IMPATOMEXPORT double get_rmsd(const core::XYZs& s0, const core::XYZs& s1);
31 
32 /**
33  Get the rmsd between two lists of particles.
34 
35  \deprecated_at{2.2} Use get_rmsd_transforming_first()
36  */
37 template <class Vector3DsOrXYZs0, class Vector3DsOrXYZs1>
38 IMPATOM_DEPRECATED_FUNCTION_DECL(2.2) inline double get_rmsd(
39  const Vector3DsOrXYZs0& m1, const Vector3DsOrXYZs1& m2,
40  const algebra::Transformation3D& tr_for_second =
41  algebra::get_identity_transformation_3d()) {
42  IMPATOM_DEPRECATED_FUNCTION_DEF(2.2, "Use it in IMP::algebra");
43  return algebra::get_rmsd_transforming_first(tr_for_second, m2, m1);
44 }
45 
46 /** RMSD on a pair of Selections.*/
47 IMPATOMEXPORT double get_rmsd(const Selection& s0, const Selection& s1);
48 
49 /** RMSD on a pair of Selections.*/
50 IMPATOMEXPORT double get_rmsd_transforming_first(
51  const algebra::Transformation3D& tr, const Selection& s0,
52  const Selection& s1);
53 
54 /**
55  \deprecated_at{2.2} Use get_rmsd_transforming_first()
56  */
57 IMPATOM_DEPRECATED_FUNCTION_DECL(2.2)
58 IMPATOMEXPORT double get_rmsd(const Selection& s0, const Selection& s1,
59  const algebra::Transformation3D& tr_for_second);
60 
61 /**
62  \deprecated_at{2.2} Use IMP::algebra::get_rmsd_transforming_first() or
63  IMP::algebra::get_rmsd()
64  */
65 IMPATOM_DEPRECATED_FUNCTION_DECL(2.2)
66 IMPATOMEXPORT double get_rmsd(const algebra::Vector3Ds& s0,
67  const algebra::Vector3Ds& s1,
68  const algebra::Transformation3D& tr_for_second);
69 
70 /**
71  \deprecated_at{2.2} Use IMP::algebra::get_rmsd_transforming_first() or
72  IMP::algebra::get_rmsd()
73  */
74 IMPATOM_DEPRECATED_FUNCTION_DECL(2.2)
75 IMPATOMEXPORT double get_rmsd(
76  const core::XYZs& s0, const core::XYZs& s1,
77  const IMP::algebra::Transformation3D& tr_for_second);
78 
79 //! Calculate the root mean square deviation between two sets of 3D points.
80 /**
81  \note the function assumes correspondence between the two sets of
82  points and does not need rigid alignment. Note that it is different from
83  get_drms().
84 
85  \genericgeometry
86  */
87 template <class Vector3DsOrXYZs0, class Vector3DsOrXYZs1>
88 inline double get_drmsd(const Vector3DsOrXYZs0& m0,
89  const Vector3DsOrXYZs1& m1) {
90  using algebra::get_vector_geometry;
91  IMP_USAGE_CHECK(m0.size() == m1.size(), "The input sets of XYZ points "
92  << "should be of the same size");
93  double drmsd = 0.0;
94 
95  int npairs = 0;
96  for (unsigned i = 0; i < m0.size() - 1; ++i) {
97  algebra::Vector3D v0i = get_vector_geometry(m0[i]);
98  algebra::Vector3D v1i = get_vector_geometry(m1[i]);
99 
100  for (unsigned j = i + 1; j < m0.size(); ++j) {
101  algebra::Vector3D v0j = get_vector_geometry(m0[j]);
102  algebra::Vector3D v1j = get_vector_geometry(m1[j]);
103 
104  double dist0 = algebra::get_distance(v0i, v0j);
105  double dist1 = algebra::get_distance(v1i, v1j);
106  drmsd += (dist0 - dist1) * (dist0 - dist1);
107  npairs++;
108  }
109  }
110  return std::sqrt(drmsd / npairs);
111 }
112 
113 //! Computes the native overlap between two sets of 3D points
114 /**
115  \param[in] m1 first set
116  \param[in] m2 second set
117  \param[in] threshold threshold distance (angstroms) for the calculation
118  \note the function assumes correspondence between two sets of points and does
119  not perform rigid alignment.
120  \return the native overlap, as a percentage (from 0 to 100)
121  \genericgeometry
122 **/
123 template <class Vector3DsOrXYZs0, class Vector3DsOrXYZs1>
124 inline double get_native_overlap(const Vector3DsOrXYZs0& m1,
125  const Vector3DsOrXYZs1& m2, double threshold) {
126  IMP_USAGE_CHECK(m1.size() == m2.size(),
127  "native_verlap: The input sets of XYZ points "
128  << "should be of the same size");
129  unsigned int distances = 0;
130  for (unsigned int i = 0; i < m1.size(); i++) {
131  double d = algebra::get_distance(get_vector_geometry(m1[i]),
132  get_vector_geometry(m2[i]));
133  if (d <= threshold) distances++;
134  }
135  return 100.0 * distances / m1.size();
136 }
137 
138 /*! Distance-RMS between two sets of points, defined as:
139  sqrt( sum[ (d1ij**2 - d2ij**2)**2]/(4 * sum[d1ij**2]) )
140  (Levitt, 1992)
141  d1ij - distance between points i and j in set m1 (the "reference" set)
142  d2ij - distance between points i and j in set m2
143  \param[in] m1 set of points
144  \param[in] m2 set of points
145 */
146 template <class Vector3DsOrXYZs0, class Vector3DsOrXYZs1>
147 inline double get_drms(const Vector3DsOrXYZs0& m1, const Vector3DsOrXYZs1& m2) {
148  IMP_USAGE_CHECK(m1.size() == m2.size(),
149  "native_verlap: The input sets of XYZ points "
150  << "should be of the same size");
151 
152  unsigned int n = m1.size();
153  double drms = 0.0;
154  double sum_d1ij = 0.0;
155  for (unsigned int i = 0; i < n; ++i) {
156  for (unsigned int j = i + 1; j < n; ++j) {
157  double sqd1 = algebra::get_squared_distance(get_vector_geometry(m1[i]),
158  get_vector_geometry(m1[j]));
159  double sqd2 = algebra::get_squared_distance(get_vector_geometry(m2[i]),
160  get_vector_geometry(m2[j]));
161  drms += (sqd1 - sqd2) * (sqd1 - sqd2);
162  sum_d1ij += sqd1;
163  }
164  }
165  drms /= (4 * sum_d1ij);
166  return sqrt(drms);
167 }
168 
169 /*! DRMS between to sets of rigid bodies. Points ij belonging to the same
170 rigid body are not evaluated, as they are the same in both sets
171  \param[in] m1 set of points
172  \param[in] m2 set of points
173  \param[in] ranges of points considered to be the same rigid body. Eg,
174  (0-29,30-49) means two rigid bodies, first with the 30 first atoms,
175  second with the last 20
176 */
177 template <class Vector3DsOrXYZs0, class Vector3DsOrXYZs1>
178 inline double get_rigid_bodies_drms(const Vector3DsOrXYZs0& m1,
179  const Vector3DsOrXYZs1& m2,
180  const IMP::IntRanges& ranges) {
181  IMP_USAGE_CHECK(m1.size() == m2.size(),
182  "get_rigid_bodies_drms: The input sets of XYZ points "
183  << "should be of the same size");
184  int n = m1.size();
185  int rn = ranges.size();
186  double drms = 0.0;
187  double sum_d1ij = 0.0;
188  for (int i = 0; i < n; ++i) {
189  int range1 = -1;
190  for (int k = 0; k < rn; ++k) {
191  if (i >= ranges[k].first && i <= ranges[k].second) {
192  range1 = k;
193  break;
194  }
195  }
196  IMP_USAGE_CHECK(range1 >= 0, "Point " << i << " of m1 does not belong to "
197  "any range");
198  for (int j = i + 1; j < n; ++j) {
199  int range2 = -1;
200  for (int k = 0; k < rn; ++k) {
201  if (j >= ranges[k].first && j <= ranges[k].second) {
202  range2 = k;
203  break;
204  }
205  }
206  IMP_USAGE_CHECK(range2 >= 0, "Point " << j << " of m2 does not belong to "
207  "any range");
208 
209  double sqd1 = algebra::get_squared_distance(get_vector_geometry(m1[i]),
210  get_vector_geometry(m1[j]));
211  sum_d1ij += sqd1;
212  if (range1 != range2) {
213  // points i and j in different ranges compare distances
214  double sqd2 = algebra::get_squared_distance(get_vector_geometry(m2[i]),
215  get_vector_geometry(m2[j]));
216  drms += (sqd1 - sqd2) * (sqd1 - sqd2);
217  }
218  }
219  }
220  drms /= (4 * sum_d1ij);
221  return sqrt(drms);
222 }
223 
224 //! Measure the difference between two placements of the same set of points
225 /**
226  \param[in] source The reference placement represented by XYZ coordinates
227  \param[in] target The modeled placement represented by XYZ coordinates
228  \note The measure quantifies the difference between placements
229  of the same structure. A rigid transformation that brings mdl1 to
230  ref1 is reported.
231  \return (d,a), A transformation from mdl to ref represented by
232  a a distance (d) and an angle (a).
233  d is the distance between the centroids of the two
234  placements and a is the axis angle of the rotation matrix between
235  the two placements
236  \note see Lasker,Topf et al JMB, 2009 for details
237  */
238 IMPATOMEXPORT FloatPair
239  get_placement_score(const core::XYZs& source, const core::XYZs& target);
240 
241 //! Measure the difference between two placements of the same set of points
242 /**
243 \param[in] ref1 The reference placement of the first component represented
244  by XYZ coordinates
245 \param[in] ref2 The reference placement of the second component represented
246  by XYZ coordinates
247 \param[in] mdl1 The modeled placement of the first component
248  represented by XYZ coordinates
249 \param[in] mdl2 The modeled placement of the second component
250  represented by XYZ coordinates
251 \return the function returns (distance score,angle score)
252 \note The measure quantifies the difference between the relative placements
253  of two components compared to a reference relative placement.
254  First, the two compared structures are brought into
255  the same frame of reference by superposing the first pair of
256  equivalent domains (ref1 and mdl1). Next, the distance and angle
257  scores are calculated for the second component using placement_score.
258 \note see Topf, Lasker et al Structure, 2008 for details
259  */
260 IMPATOMEXPORT FloatPair get_component_placement_score(const core::XYZs& ref1,
261  const core::XYZs& ref2,
262  const core::XYZs& mdl1,
263  const core::XYZs& mdl2);
264 
265 //! Measure the RMSD between two placements of the same set of points
266 /**
267 \param[in] ref1 The reference placement of the first component represented
268  by XYZ coordinates
269 \param[in] ref2 The reference placement of the second component represented
270  by XYZ coordinates
271 \param[in] mdl1 The modeled placement of the first component
272  represented by XYZ coordinates
273 \param[in] mdl2 The modeled placement of the second component
274  represented by XYZ coordinates
275 \note The measure quantifies the RMSD between the relative placements
276  of two components compared to a reference relative placement.
277  First, the two compared structures are brought into
278  the same frame of reference by superposing the first pair of
279  equivalent domains (ref1 and mdl1). Next, the RMSD is calculated
280  for the second component
281 \note see Lasker et al JMB, 2009 for details
282  */
283 IMPATOMEXPORT double get_pairwise_rmsd_score(const core::XYZs& ref1,
284  const core::XYZs& ref2,
285  const core::XYZs& mdl1,
286  const core::XYZs& mdl2);
287 
288 /** Compute the radius of gyration of a set of particles with (optional)
289  radii and mass and (non-optional) coordinates. Either all particles
290  must have mass or none of them.
291 */
292 IMPATOMEXPORT double get_radius_of_gyration(const kernel::ParticlesTemp& ps);
293 
294 IMPATOM_END_NAMESPACE
295 
296 #endif /* IMPATOM_DISTANCE_H */
double get_radius_of_gyration(const kernel::ParticlesTemp &ps)
Simple 3D transformation class.
double get_rigid_bodies_drms(const Vector3DsOrXYZs0 &m1, const Vector3DsOrXYZs1 &m2, const IMP::IntRanges &ranges)
Import IMP/kernel/base_types.h in the namespace.
double get_drms(const Vector3DsOrXYZs0 &m1, const Vector3DsOrXYZs1 &m2)
double get_native_overlap(const Vector3DsOrXYZs0 &m1, const Vector3DsOrXYZs1 &m2, double threshold)
Computes the native overlap between two sets of 3D points.
double get_drmsd(const Vector3DsOrXYZs0 &m0, const Vector3DsOrXYZs1 &m1)
Calculate the root mean square deviation between two sets of 3D points.
Definition: atom/distance.h:88
std::pair< double, double > FloatPair
A generic pair of floats.
Definition: base/types.h:27
double get_squared_distance(const VectorD< D > &v1, const VectorD< D > &v2)
compute the squared distance between two vectors
Definition: VectorD.h:201
double get_rmsd(const core::XYZs &s0, const core::XYZs &s1, const IMP::algebra::Transformation3D &tr_for_second)
Simple xyz decorator.
#define IMP_USAGE_CHECK(expr, message)
A runtime test for incorrect usage of a class or method.
Decorator for helping deal with a hierarchy of molecules.
base::Vector< VectorD< 3 > > Vector3Ds
Definition: VectorD.h:397
FloatPair get_component_placement_score(const core::XYZs &ref1, const core::XYZs &ref2, const core::XYZs &mdl1, const core::XYZs &mdl2)
Measure the difference between two placements of the same set of points.
double get_distance(const Plane3D &pln, const Vector3D &p)
Return the distance between a plane and a point in 3D.
Definition: Plane3D.h:67
Transformation3D get_identity_transformation_3d()
Return a transformation that does not do anything.
double get_pairwise_rmsd_score(const core::XYZs &ref1, const core::XYZs &ref2, const core::XYZs &mdl1, const core::XYZs &mdl2)
Measure the RMSD between two placements of the same set of points.
VectorD< 3 > Vector3D
Definition: VectorD.h:395
FloatPair get_placement_score(const core::XYZs &source, const core::XYZs &target)
Measure the difference between two placements of the same set of points.
A set of useful functionality on IMP::atom::Hierarchy decorators.
double get_rmsd_transforming_first(const algebra::Transformation3D &tr, const Selection &s0, const Selection &s1)
distance metrics