IMP logo
IMP Reference Guide  2.20.2
The Integrative Modeling Platform
rigid_fitting.h
Go to the documentation of this file.
1 /**
2  * \file IMP/em/rigid_fitting.h
3  * \brief Performs rigid fitting between a set of particles and a density map
4  *
5  * Copyright 2007-2022 IMP Inventors. All rights reserved.
6  */
7 
8 #ifndef IMPEM_RIGID_FITTING_H
9 #define IMPEM_RIGID_FITTING_H
10 
11 #include <IMP/core/MonteCarlo.h>
13 #include <IMP/RestraintSet.h>
15 #include <IMP/Particle.h>
16 #include <IMP/Model.h>
17 #include "DensityMap.h"
18 #include "FitRestraint.h"
19 #include <IMP/em/em_config.h>
20 #include <IMP/OptimizerState.h>
21 #include <IMP/core/rigid_bodies.h>
22 #include <IMP/ScoreState.h>
23 #include <algorithm>
24 #include <cereal/access.hpp>
25 #include <cereal/types/vector.hpp>
26 
27 namespace cereal {
28  template<class Archive> void serialize(Archive &ar,
29  std::pair<IMP::algebra::Transformation3D, IMP::Float> &g) {
30  ar(g.first, g.second);
31  }
32 }
33 
34 IMPEM_BEGIN_NAMESPACE
35 
36 //! A simple list of fitting solutions.
37 /**
38  \see local_rigid_fitting_around_point
39  \see local_rigid_fitting_around_points
40  \see local_rigid_fitting_grid_search
41  \see compute_fitting_scores
42  */
43 class IMPEMEXPORT FittingSolutions {
44  typedef std::pair<algebra::Transformation3D, Float> FittingSolution;
45  struct sort_by_cc {
46  bool operator()(const FittingSolution &s1,
47  const FittingSolution &s2) const {
48  return s1.second < s2.second;
49  }
50  };
51 
52  public:
53  FittingSolutions() {}
54  //! Get the number of solutions in the set
55  inline int get_number_of_solutions() const { return fs_.size(); }
56  //! Get the score of the i'th solution
57  /**
58  \return the i'th transformation, or throw an exception
59  if the index is out of range
60  */
61  inline algebra::Transformation3D get_transformation(unsigned int i) const {
62  IMP_USAGE_CHECK(i < fs_.size(),
63  "The index requested ("
64  << i << ") in get_transformation is out of range ("
65  << fs_.size() << ")" << std::endl);
66  return fs_[i].first;
67  }
68  //! Get the score of the i'th solution
69  /**
70  \return the i'th score, or throw an exception
71  if the index is out of range
72  */
73  inline Float get_score(unsigned int i) const {
74  IMP_USAGE_CHECK(i < fs_.size(),
75  "The index requested ("
76  << i << ") in get_transformation is out of range ("
77  << fs_.size() << ")" << std::endl);
78  return fs_[i].second;
79  }
80  void set_score(unsigned int i, Float score) {
81  IMP_USAGE_CHECK(i < fs_.size(),
82  "The index requested ("
83  << i << ") in get_transformation is out of range ("
84  << fs_.size() << ")" << std::endl);
85  fs_[i].second = score;
86  }
87  //! Add a solution to the fitting solution set
88  void add_solution(const algebra::Transformation3D &t, Float score);
89  //! Sort solutions by cross-correlation scores
90  void sort(bool reverse = false);
91  //! Multiply each transformation (T) by t,
92  //!such that the new transformation are T*t
94  for (unsigned int i = 0; i < fs_.size(); i++)
95  fs_[i].first = fs_[i].first * t;
96  }
97  inline algebra::Transformation3Ds get_transformations() const {
99  for (unsigned int i = 0; i < fs_.size(); i++)
100  all_ts.push_back(fs_[i].first);
101  return all_ts;
102  }
103  void show(std::ostream &out = std::cout) const {
104  for (std::vector<FittingSolution>::const_iterator it = fs_.begin();
105  it != fs_.end(); it++) {
106  out << "(" << it->first << " , " << it->second << ")" << std::endl;
107  }
108  }
109 
110  protected:
111  std::vector<FittingSolution> fs_;
112 
113 private:
114  friend class cereal::access;
115 
116  template<class Archive> void serialize(Archive &ar) {
117  ar(fs_);
118  }
119 };
121 
122 //! Local rigid fitting of a rigid body around a center point
123 /**
124 \brief Fit a set of particles to a density map around an anchor point.
125  The fitting is assessed using the cross-correlation score.
126  The optimization is a standard MC/CG procedure.
127  The function returns a list of solutions sorted by the cross-correlation
128  score.
129 \note The returned cross-correlation score is 1-cc, as we usually want to
130  minimize a scoring function. Thus a score of 1 means no-correlation
131  and a score of 0. is perfect correlation.
132 \note The input rigid body should also be IMP::atom::Hierarchy
133 \param[in] p The rigid body to fit
134 \param[in] refiner Refiner to yield rigid body members
135 \param[in] weight_key The weight key of the particles in the rigid body
136 \param[in] dmap The density map to fit to
137 \param[in] anchor_centroid The point to fit the particles around
138 \param[in] display_log If provided, then intermediate states
139  in during the optimization procedure are printed
140 \param[in] number_of_optimization_runs number of Monte Carlo optimizations
141 \param[in] number_of_mc_steps number of steps in a Monte Carlo optimization
142 \param[in] number_of_cg_steps number of Conjugate Gradients steps in
143  a Monte Carlo step
144 \param[in] max_translation maximum translation step in a MC optimization step
145 \param[in] max_rotation maximum rotation step in a single MC optimization step
146 \param[in] fast if true the density map of the rigid body is not resampled
147  but transformed at each iteration of the optimization
148 \return the refined fitting solutions
149 */
151  Particle *p, Refiner *refiner, const FloatKey &weight_key,
152  DensityMap *dmap, const algebra::Vector3D &anchor_centroid,
153  OptimizerStates display_log, Int number_of_optimization_runs = 5,
154  Int number_of_mc_steps = 10, Int number_of_cg_steps = 100,
155  Float max_translation = 2., Float max_rotation = .3, bool fast = false);
156 
157 //! Local rigid fitting of a rigid body
158 /** Fit a set of particles to a density map around their centroid.
159  The fitting is assessed using the cross-correlation score.
160  The optimization is a standard MC/CG procedure.
161  The function returns a list of solutions sorted by the cross-correlation
162  score.
163 \note The returned cross-correlation score is 1-cc, as we usually want to
164  minimize a scoring function. Thus a score of 1 means no correlation
165  and a score of 0 is perfect correlation.
166 \note The input rigid body should also be IMP::atom::Hierarchy
167 \param[in] p The root of the hierarchy to fit
168 \param[in] refiner The refiner to get the leaves of the particle
169 \param[in] weight_key The weight key of the particles in the rigid body
170 \param[in] dmap The density map to fit to
171 \param[in] display_log If provided, then intermediate states
172  during the optimization procedure are printed
173 \param[in] number_of_optimization_runs number of Monte Carlo optimizations
174 \param[in] number_of_mc_steps number of steps in a Monte Carlo optimization
175 \param[in] number_of_cg_steps number of Conjugate Gradients steps in
176  a Monte Carlo step
177 \param[in] max_translation maximum translation step in a MC optimization step
178 \param[in] max_rotation maximum rotation step in radians in a single
179  MC optimization step
180 \param[in] fast if true the density map of the rigid body is not resampled
181  but transformed at each iteration of the optimization
182 \return the refined fitting solutions
183 */
184 
186  Particle *p, Refiner *refiner, const FloatKey &weight_key,
187  DensityMap *dmap, OptimizerStates display_log,
188  Int number_of_optimization_runs = 5, Int number_of_mc_steps = 10,
189  Int number_of_cg_steps = 100, Float max_translation = 2.,
190  Float max_rotation = .3, bool fast = true) {
191  IMP_LOG_VERBOSE("Start: local_rigid_fitting\n");
192  algebra::Vector3D rb_cen =
194  IMP_LOG_VERBOSE("centroid is:" << rb_cen << "\n");
196  p, refiner, weight_key, dmap, rb_cen, display_log,
197  number_of_optimization_runs, number_of_mc_steps, number_of_cg_steps,
198  max_translation, max_rotation, fast);
199 }
200 
201 //! Local rigid fitting of a rigid body around a set of center points
202 /**
203 \brief Fit a set of particles to a density map around each of the input points.
204  The function apply local_rigid_fitting_around_point around each center.
205 \note The input rigid body should also be IMP::atom::Hierarchy
206 \param[in] p The rigid body to fit
207 \param[in] refiner Refiner to yield rigid body members
208 \param[in] wei_key The weight key of the particles in the rigid body
209 \param[in] dmap The density map to fit to
210 \param[in] anchor_centroids The points to fit the particles around
211 \param[in] display_log If provided, then intermediate states
212  in during the optimization procedure are printed
213 \param[in] number_of_optimization_runs number of Monte Carlo optimizations
214 \param[in] number_of_mc_steps number of steps in a Monte Carlo optimization
215 \param[in] number_of_cg_steps number of Conjugate Gradients steps in
216  a Monte Carlo step
217 \param[in] max_translation maximum translation step in a MC optimization step
218 \param[in] max_rotation maximum rotation step in a single MC optimization step
219 \return the refined fitting solutions
220 */
221 IMPEMEXPORT FittingSolutions local_rigid_fitting_around_points(
222  Particle *p, Refiner *refiner, const FloatKey &wei_key,
223  DensityMap *dmap, const algebra::Vector3Ds &anchor_centroids,
224  OptimizerStates display_log, Int number_of_optimization_runs = 5,
225  Int number_of_mc_steps = 10, Int number_of_cg_steps = 100,
226  Float max_translation = 2., Float max_rotation = .3);
227 
228 //! Local grid search rigid fitting
229 /**
230 \brief Fit a set of particles to a density map around their centroid.
231  The fitting is assessed using the cross-correlation score.
232  The optimization is a grid search
233 \note The transformations are not clustered.
234 \note The returned cross-correlation score is 1-cc, as we usually want to
235  minimize a scoring function. Thus a score of 1 means no-correlation
236  and a score of 0. is perfect correlation.
237 \param[in] ps The particles to be fitted (treated rigid)
238 \param[in] wei_key The weight key of the particles in the rigid body
239 \param[in] dmap The density map to fit to
240 \param[in] max_voxels_translation Sample translations within
241  -max_voxel_translation to max_voxel_translation
242 \param[in] translation_step The translation sampling step
243 \param[in] max_angle_in_radians Sample rotations with +- max_angle_in_radians
244  around the current rotation
245 \param[in] number_of_rotations The number of rotations to sample
246 \return the refined fitting solutions
247 */
248 IMPEMEXPORT FittingSolutions local_rigid_fitting_grid_search(
249  const ParticlesTemp &ps, const FloatKey &wei_key, DensityMap *dmap,
250  Int max_voxels_translation = 2, Int translation_step = 1,
251  Float max_angle_in_radians = 0.174, Int number_of_rotations = 100);
252 
253 //! Compute fitting scores for a given set of rigid transformations
254 /**
255 \brief Score how well a set of particles fit to the map in various
256  rigid transformations.
257 \param[in] ps The particles to be fitted (treated rigid)
258 \param[in] em_map The density map to fit to
259 \param[in] wei_key The weight key of the particles in the rigid body
260 \param[in] fast_version If fast_version is used the sampled density map of the
261  input particles (ps) is not resampled for each
262  transformation but instead the corresponding grid
263  is rotated. This option significantly improves the
264  running times but the returned scores are less accurate
265 \param[in] transformations A set of rigid transformations
266 \param[in] fast_version if true the density map of each transformation
267  is interpolated
268 \param[in] local_score if true a local cross correlation score is used
269 \return The scored fitting solutions
270 \note the function assumes the density map holds its density
271  */
272 IMPEMEXPORT FittingSolutions
273  compute_fitting_scores(const ParticlesTemp &ps, DensityMap *em_map,
274  const algebra::Transformation3Ds &transformations,
275  bool fast_version = false, bool local_score = false,
276  const FloatKey &wei_key =
277  atom::Mass::get_mass_key());
278 
279 //! Compute fitting scores for a given set of rigid transformations
280 /**
281 \brief Score how well a rigid body fits to the map
282 \param[in] em_map The density map to fit to
283 \param[in] rb The rigid body
284 \param[in] refiner The rigid body refiner
285 \param[in] transformations Transformations of the rigid body
286 \return The scored fitting solutions
287 \note the function assumes the density map holds its density
288  */
290  DensityMap *em_map, core::RigidBody rb, Refiner *refiner,
291  const algebra::Transformation3Ds &transformations) {
292  return compute_fitting_scores(refiner->get_refined(rb), em_map,
293  transformations, true);
294 }
295 
296 //! Compute fitting scores for a given set of rigid transformations
297 /**
298 \brief Scores how well a set of particles fit a map
299 \param[in] ps The particles to be fitted
300 \param[in] em_map The density map to fit to
301 \param[in] wei_key The weight key of the particles in the rigid body
302 \note the function assumes the density map holds its density
303  */
304 IMPEMEXPORT Float
305  compute_fitting_score(const ParticlesTemp &ps, DensityMap *em_map,
306  FloatKey wei_key = atom::Mass::get_mass_key());
307 
308 IMPEM_END_NAMESPACE
309 #endif /* IMPEM_RIGID_FITTING_H */
Simple 3D transformation class.
Simple conjugate gradients optimizer.
Calculate score based on fit to EM map.
Used to hold a set of related restraints.
algebra::Transformation3D get_transformation(unsigned int i) const
Get the score of the i'th solution.
Definition: rigid_fitting.h:61
Simple Monte Carlo optimizer.
Storage of a model, its restraints, constraints and particles.
A simple list of fitting solutions.
Definition: rigid_fitting.h:43
#define IMP_LOG_VERBOSE(expr)
Definition: log_macros.h:83
algebra::Vector3D get_centroid(const XYZs &ps)
Get the centroid.
Class for handling density maps.
Class for handling density maps.
Definition: DensityMap.h:95
#define IMP_VALUES(Name, PluralName)
Define the type for storing sets of values.
Definition: value_macros.h:23
functionality for defining rigid bodies
Float compute_fitting_score(const ParticlesTemp &ps, DensityMap *em_map, FloatKey wei_key=atom::Mass::get_mass_key())
Compute fitting scores for a given set of rigid transformations.
FittingSolutions compute_fitting_scores(DensityMap *em_map, core::RigidBody rb, Refiner *refiner, const algebra::Transformation3Ds &transformations)
Compute fitting scores for a given set of rigid transformations.
FittingSolutions local_rigid_fitting_around_point(Particle *p, Refiner *refiner, const FloatKey &weight_key, DensityMap *dmap, const algebra::Vector3D &anchor_centroid, OptimizerStates display_log, Int number_of_optimization_runs=5, Int number_of_mc_steps=10, Int number_of_cg_steps=100, Float max_translation=2., Float max_rotation=.3, bool fast=false)
Local rigid fitting of a rigid body around a center point.
Shared score state.
virtual const ParticlesTemp get_refined(Particle *a) const =0
Refine the passed particle into a set of particles.
Simple 3D transformation class.
Classes to handle individual model particles. (Note that implementation of inline functions is in int...
std::ostream & show(Hierarchy h, std::ostream &out=std::cout)
Print the hierarchy using a given decorator to display each node.
FittingSolutions local_rigid_fitting(Particle *p, Refiner *refiner, const FloatKey &weight_key, DensityMap *dmap, OptimizerStates display_log, Int number_of_optimization_runs=5, Int number_of_mc_steps=10, Int number_of_cg_steps=100, Float max_translation=2., Float max_rotation=.3, bool fast=true)
Local rigid fitting of a rigid body.
VectorD< 3 > Vector3D
Definition: VectorD.h:425
double Float
Basic floating-point value (could be float, double...)
Definition: types.h:19
Abstract class to implement hierarchical methods.
Definition: Refiner.h:34
void multiply(const algebra::Transformation3D &t)
Definition: rigid_fitting.h:93
Class to handle individual particles of a Model object.
Definition: Particle.h:43
#define IMP_USAGE_CHECK(expr, message)
A runtime test for incorrect usage of a class or method.
Definition: check_macros.h:168
int Int
Basic integer value.
Definition: types.h:34
Shared optimizer state.
FittingSolutions local_rigid_fitting_around_points(Particle *p, Refiner *refiner, const FloatKey &wei_key, DensityMap *dmap, const algebra::Vector3Ds &anchor_centroids, OptimizerStates display_log, Int number_of_optimization_runs=5, Int number_of_mc_steps=10, Int number_of_cg_steps=100, Float max_translation=2., Float max_rotation=.3)
Local rigid fitting of a rigid body around a set of center points.
FittingSolutions local_rigid_fitting_grid_search(const ParticlesTemp &ps, const FloatKey &wei_key, DensityMap *dmap, Int max_voxels_translation=2, Int translation_step=1, Float max_angle_in_radians=0.174, Int number_of_rotations=100)
Local grid search rigid fitting.
A decorator for a rigid body.
Definition: rigid_bodies.h:82
int get_number_of_solutions() const
Get the number of solutions in the set.
Definition: rigid_fitting.h:55
Float get_score(unsigned int i) const
Get the score of the i'th solution.
Definition: rigid_fitting.h:73