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