IMP  2.4.0
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-2015 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 distance 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 
114 //! Calculate the distance root mean square deviation between two sets of 3D points
115 //! with a distance threshold
116 /**
117  \note the function assumes correspondence between the two sets of
118  points and does not need rigid alignment. Note that it is different from
119  get_drmsd().
120 
121  \genericgeometry
122  */
123 template <class Vector3DsOrXYZs0, class Vector3DsOrXYZs1>
124 inline double get_drmsd_Q(const Vector3DsOrXYZs0& m0,
125  const Vector3DsOrXYZs1& m1, double threshold) {
126  using algebra::get_vector_geometry;
127  IMP_USAGE_CHECK(m0.size() == m1.size(), "The input sets of XYZ points "
128  << "should be of the same size");
129  double drmsd = 0.0;
130 
131  int npairs = 0;
132  for (unsigned i = 0; i < m0.size() - 1; ++i) {
133  algebra::Vector3D v0i = get_vector_geometry(m0[i]);
134  algebra::Vector3D v1i = get_vector_geometry(m1[i]);
135 
136  for (unsigned j = i + 1; j < m0.size(); ++j) {
137  algebra::Vector3D v0j = get_vector_geometry(m0[j]);
138  algebra::Vector3D v1j = get_vector_geometry(m1[j]);
139 
140  double dist0 = algebra::get_distance(v0i, v0j);
141  double dist1 = algebra::get_distance(v1i, v1j);
142  if (dist0<=threshold || dist1<=threshold){
143  drmsd += (dist0 - dist1) * (dist0 - dist1);
144  npairs++;
145  }
146  }
147  }
148  return std::sqrt(drmsd / npairs);
149 }
150 
151 
152 //! Computes the native overlap between two sets of 3D points
153 /**
154  \param[in] m1 first set
155  \param[in] m2 second set
156  \param[in] threshold threshold distance (angstroms) for the calculation
157  \note the function assumes correspondence between two sets of points and does
158  not perform rigid alignment.
159  \return the native overlap, as a percentage (from 0 to 100)
160  \genericgeometry
161 **/
162 template <class Vector3DsOrXYZs0, class Vector3DsOrXYZs1>
163 inline double get_native_overlap(const Vector3DsOrXYZs0& m1,
164  const Vector3DsOrXYZs1& m2, double threshold) {
165  IMP_USAGE_CHECK(m1.size() == m2.size(),
166  "native_overlap: The input sets of XYZ points "
167  << "should be of the same size");
168  unsigned int distances = 0;
169  for (unsigned int i = 0; i < m1.size(); i++) {
170  double d = algebra::get_distance(get_vector_geometry(m1[i]),
171  get_vector_geometry(m2[i]));
172  if (d <= threshold) distances++;
173  }
174  return 100.0 * distances / m1.size();
175 }
176 
177 /*! Distance-RMS between two sets of points, defined as:
178  sqrt( sum[ (d1ij**2 - d2ij**2)**2]/(4 * sum[d1ij**2]) )
179  (Levitt, 1992)
180  d1ij - distance between points i and j in set m1 (the "reference" set)
181  d2ij - distance between points i and j in set m2
182  \param[in] m1 set of points
183  \param[in] m2 set of points
184 */
185 template <class Vector3DsOrXYZs0, class Vector3DsOrXYZs1>
186 inline double get_drms(const Vector3DsOrXYZs0& m1, const Vector3DsOrXYZs1& m2) {
187  IMP_USAGE_CHECK(m1.size() == m2.size(),
188  "native_overlap: The input sets of XYZ points "
189  << "should be of the same size");
190 
191  unsigned int n = m1.size();
192  double drms = 0.0;
193  double sum_d1ij = 0.0;
194  for (unsigned int i = 0; i < n; ++i) {
195  for (unsigned int j = i + 1; j < n; ++j) {
196  double sqd1 = algebra::get_squared_distance(get_vector_geometry(m1[i]),
197  get_vector_geometry(m1[j]));
198  double sqd2 = algebra::get_squared_distance(get_vector_geometry(m2[i]),
199  get_vector_geometry(m2[j]));
200  drms += (sqd1 - sqd2) * (sqd1 - sqd2);
201  sum_d1ij += sqd1;
202  }
203  }
204  drms /= (4 * sum_d1ij);
205  return sqrt(drms);
206 }
207 
208 /*! DRMS between to sets of rigid bodies. Points ij belonging to the same
209 rigid body are not evaluated, as they are the same in both sets
210  \param[in] m1 set of points
211  \param[in] m2 set of points
212  \param[in] ranges of points considered to be the same rigid body. Eg,
213  (0-29,30-49) means two rigid bodies, first with the 30 first atoms,
214  second with the last 20
215 */
216 template <class Vector3DsOrXYZs0, class Vector3DsOrXYZs1>
217 inline double get_rigid_bodies_drms(const Vector3DsOrXYZs0& m1,
218  const Vector3DsOrXYZs1& m2,
219  const IMP::IntRanges& ranges) {
220  IMP_USAGE_CHECK(m1.size() == m2.size(),
221  "get_rigid_bodies_drms: The input sets of XYZ points "
222  << "should be of the same size");
223  int n = m1.size();
224  int rn = ranges.size();
225  double drms = 0.0;
226  double sum_d1ij = 0.0;
227  for (int i = 0; i < n; ++i) {
228  int range1 = -1;
229  for (int k = 0; k < rn; ++k) {
230  if (i >= ranges[k].first && i <= ranges[k].second) {
231  range1 = k;
232  break;
233  }
234  }
235  IMP_USAGE_CHECK(range1 >= 0, "Point " << i << " of m1 does not belong to "
236  "any range");
237  for (int j = i + 1; j < n; ++j) {
238  int range2 = -1;
239  for (int k = 0; k < rn; ++k) {
240  if (j >= ranges[k].first && j <= ranges[k].second) {
241  range2 = k;
242  break;
243  }
244  }
245  IMP_USAGE_CHECK(range2 >= 0, "Point " << j << " of m2 does not belong to "
246  "any range");
247 
248  double sqd1 = algebra::get_squared_distance(get_vector_geometry(m1[i]),
249  get_vector_geometry(m1[j]));
250  sum_d1ij += sqd1;
251  if (range1 != range2) {
252  // points i and j in different ranges compare distances
253  double sqd2 = algebra::get_squared_distance(get_vector_geometry(m2[i]),
254  get_vector_geometry(m2[j]));
255  drms += (sqd1 - sqd2) * (sqd1 - sqd2);
256  }
257  }
258  }
259  drms /= (4 * sum_d1ij);
260  return sqrt(drms);
261 }
262 
263 //! Measure the difference between two placements of the same set of points
264 /**
265  \param[in] source The reference placement represented by XYZ coordinates
266  \param[in] target The modeled placement represented by XYZ coordinates
267  \note The measure quantifies the difference between placements
268  of the same structure. A rigid transformation that brings mdl1 to
269  ref1 is reported.
270  \return (d,a), A transformation from mdl to ref represented by
271  a a distance (d) and an angle (a).
272  d is the distance between the centroids of the two
273  placements and a is the axis angle of the rotation matrix between
274  the two placements
275  \note see Lasker,Topf et al JMB, 2009 for details
276  */
277 IMPATOMEXPORT FloatPair
278  get_placement_score(const core::XYZs& source, const core::XYZs& target);
279 
280 //! Measure the difference between two placements of the same set of points
281 /**
282 \param[in] ref1 The reference placement of the first component represented
283  by XYZ coordinates
284 \param[in] ref2 The reference placement of the second component represented
285  by XYZ coordinates
286 \param[in] mdl1 The modeled placement of the first component
287  represented by XYZ coordinates
288 \param[in] mdl2 The modeled placement of the second component
289  represented by XYZ coordinates
290 \return the function returns (distance score,angle score)
291 \note The measure quantifies the difference between the relative placements
292  of two components compared to a reference relative placement.
293  First, the two compared structures are brought into
294  the same frame of reference by superposing the first pair of
295  equivalent domains (ref1 and mdl1). Next, the distance and angle
296  scores are calculated for the second component using placement_score.
297 \note see Topf, Lasker et al Structure, 2008 for details
298  */
299 IMPATOMEXPORT FloatPair get_component_placement_score(const core::XYZs& ref1,
300  const core::XYZs& ref2,
301  const core::XYZs& mdl1,
302  const core::XYZs& mdl2);
303 
304 //! Measure the RMSD between two placements of the same set of points
305 /**
306 \param[in] ref1 The reference placement of the first component represented
307  by XYZ coordinates
308 \param[in] ref2 The reference placement of the second component represented
309  by XYZ coordinates
310 \param[in] mdl1 The modeled placement of the first component
311  represented by XYZ coordinates
312 \param[in] mdl2 The modeled placement of the second component
313  represented by XYZ coordinates
314 \note The measure quantifies the RMSD between the relative placements
315  of two components compared to a reference relative placement.
316  First, the two compared structures are brought into
317  the same frame of reference by superposing the first pair of
318  equivalent domains (ref1 and mdl1). Next, the RMSD is calculated
319  for the second component
320 \note see Lasker et al JMB, 2009 for details
321  */
322 IMPATOMEXPORT double get_pairwise_rmsd_score(const core::XYZs& ref1,
323  const core::XYZs& ref2,
324  const core::XYZs& mdl1,
325  const core::XYZs& mdl2);
326 
327 /** Compute the radius of gyration of a set of particles with (optional)
328  radii and mass and (non-optional) coordinates. Either all particles
329  must have mass or none of them.
330 */
331 IMPATOMEXPORT double get_radius_of_gyration(const kernel::ParticlesTemp& ps);
332 
333 IMPATOM_END_NAMESPACE
334 
335 #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 distance 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: 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.
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:71
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
#define IMP_USAGE_CHECK(expr, message)
A runtime test for incorrect usage of a class or method.
Definition: check_macros.h:170
FloatPair get_placement_score(const core::XYZs &source, const core::XYZs &target)
Measure the difference between two placements of the same set of points.
double get_drmsd_Q(const Vector3DsOrXYZs0 &m0, const Vector3DsOrXYZs1 &m1, double threshold)
Select a subset of a hierarchy.
double get_rmsd_transforming_first(const algebra::Transformation3D &tr, const Selection &s0, const Selection &s1)
distance metrics