IMP logo
IMP Reference Guide  2.14.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-2020 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  (currently unused)
52  \param[in] model_cutoff_dist Cutoff for the model-model interactions
53  \param[in] density_cutoff_dist Cutoff for model-density interactions
54  \param[in] slope Gentle term to move all particles to the density
55  \param[in] update_model (DEPRECATED) update model each cycle
56  \param[in] backbone_slope Limit the slope only to backbone particles
57  \param[in] local Only consider density particles that are within the
58  specified cutoff of the model particles (experimental)
59  \param[in] name Name of this restraint
60  \note the model and density particles must be set up as Gaussian
61  */
63  ParticleIndexes model_ps, ParticleIndexes density_ps,
64  ParticleIndex global_sigma,
65  Float model_cutoff_dist, Float density_cutoff_dist,
66  Float slope,
67  bool update_model=true, bool backbone_slope=false,
68  bool local=false,
69  std::string name="GaussianEMRestraint%1%");
70 
71  //! Returns exp(score)
72  double get_probability() const {
73  return exp(-unprotected_evaluate(NULL));
74  }
75 
76  //! Get cross correlation between the model and the map
77  /** This CCC is that calculated from the last scoring function
78  evaluation; calling this function before the score is calculated
79  results in undefined behavior.
80  */
82  return cross_correlation_;
83  }
84 
85  //! Set the filename corresponding to the density GMM particles
86  /** If the density GMM particles were read from a file, this method
87  can be used to tell the restraint so that it can track this
88  information back to the original EM density file, which is useful
89  for deposition. */
90  void set_density_filename(std::string density_fn) {
91  density_fn_ = get_absolute_path(density_fn);
92  }
93 
94  //! Pre-calculate the density-density and model-model scores
95  /** This is automatically called by the constructor.
96  You only need to call it manually if you change Gaussian variances.
97  */
98  void compute_initial_scores();
99 
100  //! Set restraint slope
101  void set_slope(Float s){slope_=s;}
102 
103  //! Get restraint slope
104  Float get_slope(){return slope_;}
105  virtual double
106  unprotected_evaluate(IMP::DerivativeAccumulator *accum) const IMP_OVERRIDE;
108  void show(std::ostream &out) const { out << "GEM restraint"; }
109 
110  //! \return Information for writing to RMF files
112 
113  //! \return Information for writing to RMF files
114  RestraintInfo *get_dynamic_info() const IMP_OVERRIDE;
115 
117 
118  private:
119  double model_cutoff_dist_, density_cutoff_dist_;
120  ParticleIndexes model_ps_;
121  ParticleIndexes density_ps_;
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:35
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