IMP logo
IMP Reference Guide  develop.cb469e0f11,2024/05/28
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-2022 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 /** RMSD on a pair of Selections.*/
33 IMPATOMEXPORT double get_rmsd(const Selection& s0, const Selection& s1);
34 
35 /** RMSD on a pair of Selections.*/
36 IMPATOMEXPORT double get_rmsd_transforming_first(
37  const algebra::Transformation3D& tr, const Selection& s0,
38  const Selection& s1);
39 
40 //! Calculate distance the root mean square deviation between two sets of 3D points.
41 /**
42  \note the function assumes correspondence between the two sets of
43  points and does not need rigid alignment. Note that it is different from
44  get_drms().
45 
46  \genericgeometry
47  */
48 template <class Vector3DsOrXYZs0, class Vector3DsOrXYZs1>
49 inline double get_drmsd(const Vector3DsOrXYZs0& m0,
50  const Vector3DsOrXYZs1& m1) {
51  using algebra::get_vector_geometry;
52  IMP_USAGE_CHECK(m0.size() == m1.size(), "The input sets of XYZ points "
53  << "should be of the same size");
54  double drmsd = 0.0;
55 
56  int npairs = 0;
57  for (unsigned i = 0; i < m0.size() - 1; ++i) {
58  algebra::Vector3D v0i = get_vector_geometry(m0[i]);
59  algebra::Vector3D v1i = get_vector_geometry(m1[i]);
60 
61  for (unsigned j = i + 1; j < m0.size(); ++j) {
62  algebra::Vector3D v0j = get_vector_geometry(m0[j]);
63  algebra::Vector3D v1j = get_vector_geometry(m1[j]);
64 
65  double dist0 = algebra::get_distance(v0i, v0j);
66  double dist1 = algebra::get_distance(v1i, v1j);
67  drmsd += (dist0 - dist1) * (dist0 - dist1);
68  npairs++;
69  }
70  }
71  return std::sqrt(drmsd / npairs);
72 }
73 
74 
75 //! Calculate the distance root mean square deviation between two sets of 3D points
76 //! with a distance threshold
77 /**
78  \note the function assumes correspondence between the two sets of
79  points and does not need rigid alignment. Note that it is different from
80  get_drmsd().
81 
82  \genericgeometry
83  */
84 template <class Vector3DsOrXYZs0, class Vector3DsOrXYZs1>
85 inline double get_drmsd_Q(const Vector3DsOrXYZs0& m0,
86  const Vector3DsOrXYZs1& m1, double threshold) {
87  using algebra::get_vector_geometry;
88  IMP_USAGE_CHECK(m0.size() == m1.size(), "The input sets of XYZ points "
89  << "should be of the same size");
90  double drmsd = 0.0;
91 
92  int npairs = 0;
93  for (unsigned i = 0; i < m0.size() - 1; ++i) {
94  algebra::Vector3D v0i = get_vector_geometry(m0[i]);
95  algebra::Vector3D v1i = get_vector_geometry(m1[i]);
96 
97  for (unsigned j = i + 1; j < m0.size(); ++j) {
98  algebra::Vector3D v0j = get_vector_geometry(m0[j]);
99  algebra::Vector3D v1j = get_vector_geometry(m1[j]);
100 
101  double dist0 = algebra::get_distance(v0i, v0j);
102  double dist1 = algebra::get_distance(v1i, v1j);
103  if (dist0<=threshold || dist1<=threshold){
104  drmsd += (dist0 - dist1) * (dist0 - dist1);
105  npairs++;
106  }
107  }
108  }
109  return std::sqrt(drmsd / npairs);
110 }
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_overlap: 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_overlap: 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 two sets of rigid bodies.
170 /** Points ij belonging to the same rigid body are not evaluated, as
171  they are the same in both sets
172  \param[in] m1 set of points
173  \param[in] m2 set of points
174  \param[in] ranges of points considered to be the same rigid body. Eg,
175  (0-29,30-49) means two rigid bodies, first with the 30 first atoms,
176  second with the last 20
177 */
178 template <class Vector3DsOrXYZs0, class Vector3DsOrXYZs1>
179 inline double get_rigid_bodies_drms(const Vector3DsOrXYZs0& m1,
180  const Vector3DsOrXYZs1& m2,
181  const IMP::IntRanges& ranges) {
182  IMP_USAGE_CHECK(m1.size() == m2.size(),
183  "get_rigid_bodies_drms: The input sets of XYZ points "
184  << "should be of the same size");
185  int n = m1.size();
186  int rn = ranges.size();
187  double drms = 0.0;
188  double sum_d1ij = 0.0;
189  for (int i = 0; i < n; ++i) {
190  int range1 = -1;
191  for (int k = 0; k < rn; ++k) {
192  if (i >= ranges[k].first && i <= ranges[k].second) {
193  range1 = k;
194  break;
195  }
196  }
197  IMP_USAGE_CHECK(range1 >= 0, "Point " << i << " of m1 does not belong to "
198  "any range");
199  for (int j = i + 1; j < n; ++j) {
200  int range2 = -1;
201  for (int k = 0; k < rn; ++k) {
202  if (j >= ranges[k].first && j <= ranges[k].second) {
203  range2 = k;
204  break;
205  }
206  }
207  IMP_USAGE_CHECK(range2 >= 0, "Point " << j << " of m2 does not belong to "
208  "any range");
209 
210  double sqd1 = algebra::get_squared_distance(get_vector_geometry(m1[i]),
211  get_vector_geometry(m1[j]));
212  sum_d1ij += sqd1;
213  if (range1 != range2) {
214  // points i and j in different ranges compare distances
215  double sqd2 = algebra::get_squared_distance(get_vector_geometry(m2[i]),
216  get_vector_geometry(m2[j]));
217  drms += (sqd1 - sqd2) * (sqd1 - sqd2);
218  }
219  }
220  }
221  drms /= (4 * sum_d1ij);
222  return sqrt(drms);
223 }
224 
225 //! Measure the difference between two placements of the same set of points
226 /**
227  \param[in] source The reference placement represented by XYZ coordinates
228  \param[in] target The modeled placement represented by XYZ coordinates
229  \note The measure quantifies the difference between placements
230  of the same structure. A rigid transformation that brings mdl1 to
231  ref1 is reported.
232  \return (d,a), A transformation from mdl to ref represented by
233  a a distance (d) and an angle (a).
234  d is the distance between the centroids of the two
235  placements and a is the axis angle of the rotation matrix between
236  the two placements
237  \note see Lasker,Topf et al JMB, 2009 for details
238  */
239 IMPATOMEXPORT FloatPair
240  get_placement_score(const core::XYZs& source, const core::XYZs& target);
241 
242 //! Measure the difference between two placements of the same set of points
243 /**
244 \param[in] ref1 The reference placement of the first component represented
245  by XYZ coordinates
246 \param[in] ref2 The reference placement of the second component represented
247  by XYZ coordinates
248 \param[in] mdl1 The modeled placement of the first component
249  represented by XYZ coordinates
250 \param[in] mdl2 The modeled placement of the second component
251  represented by XYZ coordinates
252 \return the function returns (distance score,angle score)
253 \note The measure quantifies the difference between the relative placements
254  of two components compared to a reference relative placement.
255  First, the two compared structures are brought into
256  the same frame of reference by superposing the first pair of
257  equivalent domains (ref1 and mdl1). Next, the distance and angle
258  scores are calculated for the second component using placement_score.
259 \note see Topf, Lasker et al Structure, 2008 for details
260  */
261 IMPATOMEXPORT FloatPair get_component_placement_score(const core::XYZs& ref1,
262  const core::XYZs& ref2,
263  const core::XYZs& mdl1,
264  const core::XYZs& mdl2);
265 
266 //! Measure the RMSD between two placements of the same set of points
267 /**
268 \param[in] ref1 The reference placement of the first component represented
269  by XYZ coordinates
270 \param[in] ref2 The reference placement of the second component represented
271  by XYZ coordinates
272 \param[in] mdl1 The modeled placement of the first component
273  represented by XYZ coordinates
274 \param[in] mdl2 The modeled placement of the second component
275  represented by XYZ coordinates
276 \note The measure quantifies the RMSD between the relative placements
277  of two components compared to a reference relative placement.
278  First, the two compared structures are brought into
279  the same frame of reference by superposing the first pair of
280  equivalent domains (ref1 and mdl1). Next, the RMSD is calculated
281  for the second component
282 \note see Lasker et al JMB, 2009 for details
283  */
284 IMPATOMEXPORT
286 (const core::XYZs& ref1,
287  const core::XYZs& ref2,
288  const core::XYZs& mdl1,
289  const core::XYZs& mdl2);
290 
291 /** Compute the radius of gyration $R_g$ of a set of particles. It is assumed
292  the particles has coordinates (XYZ decorator), and optionally, the $R_g$
293  may be weighted by mass or radius^3.
294 
295  @param ps the particles
296  @param weighted If true then $R_g$ is weighted by either mass or radius^3,
297  in this order of preference. If none exists, an unweighted $R_g$ is
298  computed. If weighted is true, then it is is assumed that either
299  all particles have mass/radius or none of them do.
300 
301  @return the weighted or unweighted radius of gyration of ps
302 */
303 IMPATOMEXPORT
305 (const ParticlesTemp& ps,
306  bool weighted=true);
307 
308 IMPATOM_END_NAMESPACE
309 
310 #endif /* IMPATOM_DISTANCE_H */
Simple 3D transformation class.
double get_rigid_bodies_drms(const Vector3DsOrXYZs0 &m1, const Vector3DsOrXYZs1 &m2, const IMP::IntRanges &ranges)
DRMS between two sets of rigid bodies.
Basic types used by IMP.
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:49
std::pair< double, double > FloatPair
A generic pair of floats.
Definition: types.h:26
double get_squared_distance(const VectorD< D > &v1, const VectorD< D > &v2)
Compute the squared distance between two vectors.
Definition: VectorD.h:188
A more IMP-like version of the std::vector.
Definition: Vector.h:50
Simple XYZ decorator.
Decorator for helping deal with a hierarchy of molecules.
double get_rmsd(const Selection &s0, const Selection &s1)
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_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:408
#define IMP_USAGE_CHECK(expr, message)
A runtime test for incorrect usage of a class or method.
Definition: check_macros.h:168
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_radius_of_gyration(const ParticlesTemp &ps, bool weighted=true)
double get_drmsd_Q(const Vector3DsOrXYZs0 &m0, const Vector3DsOrXYZs1 &m1, double threshold)
Definition: atom/distance.h:85
Select a subset of a hierarchy.
double get_rmsd_transforming_first(const algebra::Transformation3D &tr, const Selection &s0, const Selection &s1)
Distance metrics.
double get_distance(const Line3D &s, const Vector3D &p)
Get closest distance between a line and a point.