IMP logo
IMP Reference Guide  2.12.0
The Integrative Modeling Platform
isd/GaussianEMRestraint.h
Go to the documentation of this file.
1 /**
2  * \file IMP/isd/GaussianEMRestraint.h
3  * \brief Restrain two sets of Gaussians (model and GMM derived from EM map)
4  *
5  * Copyright 2007-2019 IMP Inventors. All rights reserved.
6  *
7  */
8 
9 #ifndef IMPISD_GAUSSIAN_EM_RESTRAINT_H
10 #define IMPISD_GAUSSIAN_EM_RESTRAINT_H
11 
12 #include "isd_config.h"
13 #include <IMP/file.h>
14 #include <IMP/PairContainer.h>
16 #include <IMP/container_macros.h>
17 #include <IMP/core/XYZ.h>
18 #include <IMP/core/Gaussian.h>
19 #include <IMP/algebra/Gaussian3D.h>
22 #include <IMP/atom/Mass.h>
23 #include <math.h>
24 #include <Eigen/Dense>
25 #include <boost/unordered_map.hpp>
26 
27 IMPISD_BEGIN_NAMESPACE
28 
29 //! Restraint between two Gaussian Mixture Models, "model" and "density"
30 /** This restrains two sets of GMMs with a density overlap function.
31  The function is correlation of the two GMMs \f$f_M\f$ and \f$f_D\f$:
32  \f[
33  \frac{2\int{f_M(x)f_D(x)dx}}{\int{f_M^2(x)+f_D^2(x)dx}}
34  \f]
35  Where the integral is the "overlap function" given by:
36  \f[
37  ov(f_M,f_D) = \sum_{i=1}^{N_M} \sum_{j=1}^{N_D} \frac{1}{(2 \pi)^{3/2}|\Sigma_{Mi}+\Sigma_{Dj}|^{1/2}}\exp\left [-\frac{1}{2}(\boldsymbol\mu_{Mi} - \boldsymbol\mu_{Dj})^\top (\Sigma_{Mi}+\Sigma_{Dj})^{-1} (\boldsymbol\mu_{Mi} - \boldsymbol \mu_{Dj})\right ]
38  \f]
39  \note Source: Greenberg, Pellarin, Sali. In preparation.
40  */
41 class IMPISDEXPORT GaussianEMRestraint : public Restraint
42 {
43 
44  public:
45  //! Setup the GaussianEMRestraint
46  /**
47  \param[in] mdl the Model object to operate on
48  \param[in] model_ps particles for the model GMM
49  \param[in] density_ps particles for the density GMM
50  \param[in] global_sigma Particle to modulate the uncertainty
51  \param[in] model_cutoff_dist Cutoff for the model-model interactions
52  \param[in] density_cutoff_dist Cutoff for model-density interactions
53  \param[in] slope Gentle term to move all particles to the density
54  \param[in] update_model (DEPRECATED) update model each cycle
55  \param[in] backbone_slope Limit the slope only to backbone particles
56  \param[in] local Only consider density particles that are within the
57  specified cutoff of the model particles (experimental)
58  \param[in] name Name of this restraint
59  \note the model and density particles must be set up as Gaussian
60  */
62  ParticleIndexes model_ps, ParticleIndexes density_ps,
63  ParticleIndex global_sigma,
64  Float model_cutoff_dist, Float density_cutoff_dist,
65  Float slope,
66  bool update_model=true, bool backbone_slope=false,
67  bool local=false,
68  std::string name="GaussianEMRestraint%1%");
69 
70  //! Returns exp(score)
71  double get_probability() const {
72  return exp(-unprotected_evaluate(NULL));
73  }
74 
75  //! Get cross correlation between the model and the map
76  /** This CCC is that calculated from the last scoring function
77  evaluation; calling this function before the score is calculated
78  results in undefined behavior.
79  */
81  return cross_correlation_;
82  }
83 
84  //! Set the filename corresponding to the density GMM particles
85  /** If the density GMM particles were read from a file, this method
86  can be used to tell the restraint so that it can track this
87  information back to the original EM density file, which is useful
88  for deposition. */
89  void set_density_filename(std::string density_fn) {
90  density_fn_ = get_absolute_path(density_fn);
91  }
92 
93  //! Pre-calculate the density-density and model-model scores
94  /** This is automatically called by the constructor.
95  You only need to call it manually if you change Gaussian variances.
96  */
97  void compute_initial_scores();
98 
99  //! Set restraint slope
100  void set_slope(Float s){slope_=s;}
101 
102  //! Get restraint slope
103  Float get_slope(){return slope_;}
104  virtual double
105  unprotected_evaluate(IMP::DerivativeAccumulator *accum) const IMP_OVERRIDE;
107  void show(std::ostream &out) const { out << "GEM restraint"; }
108 
109  //! \return Information for writing to RMF files
111 
112  //! \return Information for writing to RMF files
113  RestraintInfo *get_dynamic_info() const IMP_OVERRIDE;
114 
116 
117  private:
118  double model_cutoff_dist_, density_cutoff_dist_;
119  ParticleIndexes model_ps_;
120  ParticleIndexes density_ps_;
121  ParticleIndex global_sigma_;
122  Float slope_;
123  bool update_model_,local_;
124  int msize_,dsize_;
125  Float normalization_;
126  Float dd_score_;
127  Float self_mm_score_;
128  PointerMember<container::CloseBipartitePairContainer> md_container_;
129  PointerMember<container::ClosePairContainer> mm_container_;
130  ParticleIndexes slope_ps_; //experiment
131  std::string density_fn_;
132 
133  //variables needed to tabulate the exponential
134  Floats exp_grid_;
135  double invdx_;
136  double argmax_;
137 
138  mutable double cross_correlation_;
139 };
140 
141 IMPISD_END_NAMESPACE
142 
143 #endif /* IMPISD_GAUSSIAN_EM_RESTRAINT_H */
Declare an efficient stl-compatible map.
Decorator to hold Gaussian3D.
A decorator for particles with mass.
virtual RestraintInfo * get_static_info() const
Definition: Restraint.h:117
Store a list of ParticleIndexes.
Return all pairs from a SingletonContainer.
#define IMP_OBJECT_METHODS(Name)
Define the basic things needed by any Object.
Definition: object_macros.h:25
Float get_slope()
Get restraint slope.
void set_density_filename(std::string density_fn)
Set the filename corresponding to the density GMM particles.
Handling of file input/output.
void set_slope(Float s)
Set restraint slope.
Macros to define containers of objects.
Simple XYZ decorator.
Class for storing model, its restraints, constraints, and particles.
Definition: Model.h:72
double get_probability() const
Returns exp(score)
A container for Pairs.
double get_cross_correlation_coefficient() const
Get cross correlation between the model and the map.
Restraint between two Gaussian Mixture Models, "model" and "density".
Return all pairs from a SingletonContainer.
A smart pointer to a ref-counted Object that is a class member.
Definition: Pointer.h:146
Report key:value information on restraints.
Definition: RestraintInfo.h:24
std::ostream & show(Hierarchy h, std::ostream &out=std::cout)
Print the hierarchy using a given decorator to display each node.
std::string get_absolute_path(std::string file)
Convert a possibly relative path to an absolute path.
double Float
Basic floating-point value (could be float, double...)
Definition: types.h:20
Gaussian shape.
virtual ModelObjectsTemp do_get_inputs() const =0
#define IMP_OVERRIDE
Cause a compile error if this method does not override a parent method.
Class for adding derivatives from restraints to the model.
A restraint is a term in an IMP ScoringFunction.
Definition: Restraint.h:54