IMP logo
IMP Reference Guide  develop.374319d7bb,2024/06/17
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-2022 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 <IMP/isd/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 <cmath>
24 #include <Eigen/Dense>
25 #include <boost/unordered_map.hpp>
26 #include <cereal/access.hpp>
27 #include <cereal/types/vector.hpp>
28 
29 IMPISD_BEGIN_NAMESPACE
30 
31 //! Restraint between two Gaussian Mixture Models, "model" and "density"
32 /** This restrains two sets of GMMs with a density overlap function.
33  The function is correlation of the two GMMs \f$f_M\f$ and \f$f_D\f$:
34  \f[
35  \frac{2\int{f_M(x)f_D(x)dx}}{\int{f_M^2(x)+f_D^2(x)dx}}
36  \f]
37  Where the integral is the "overlap function" given by:
38  \f[
39  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 ]
40  \f]
41  \note Source: Greenberg, Pellarin, Sali. In preparation.
42  */
43 class IMPISDEXPORT GaussianEMRestraint : public Restraint
44 {
45 
46  public:
47  //! Setup the GaussianEMRestraint
48  /**
49  \param[in] mdl the Model object to operate on
50  \param[in] model_ps particles for the model GMM
51  \param[in] density_ps particles for the density GMM
52  \param[in] global_sigma Particle to modulate the uncertainty
53  (currently unused)
54  \param[in] model_cutoff_dist Cutoff for the model-model interactions
55  \param[in] density_cutoff_dist Cutoff for model-density interactions
56  \param[in] slope Gentle term to move all particles to the density
57  \param[in] update_model (DEPRECATED) update model each cycle
58  \param[in] backbone_slope Limit the slope only to backbone particles
59  \param[in] local Only consider density particles that are within the
60  specified cutoff of the model particles (experimental)
61  \param[in] name Name of this restraint
62  \note the model and density particles must be set up as Gaussian
63  */
65  ParticleIndexes model_ps, ParticleIndexes density_ps,
66  ParticleIndex global_sigma,
67  Float model_cutoff_dist, Float density_cutoff_dist,
68  Float slope,
69  bool update_model=true, bool backbone_slope=false,
70  bool local=false,
71  std::string name="GaussianEMRestraint%1%");
73 
74  //! Returns exp(score)
75  double get_probability() const {
76  return exp(-unprotected_evaluate(NULL));
77  }
78 
79  //! Get cross correlation between the model and the map
80  /** This CCC is that calculated from the last scoring function
81  evaluation; calling this function before the score is calculated
82  results in undefined behavior.
83  */
85  return cross_correlation_;
86  }
87 
88  //! Set the filename corresponding to the density GMM particles
89  /** If the density GMM particles were read from a file, this method
90  can be used to tell the restraint so that it can track this
91  information back to the original EM density file, which is useful
92  for deposition. */
93  void set_density_filename(std::string density_fn) {
94  density_fn_ = get_absolute_path(density_fn);
95  }
96 
97  //! Pre-calculate the density-density and model-model scores
98  /** This is automatically called by the constructor.
99  You only need to call it manually if you change Gaussian variances.
100  */
101  void compute_initial_scores();
102 
103  //! Set restraint slope
104  void set_slope(Float s){slope_=s;}
105 
106  //! Get restraint slope
107  Float get_slope(){return slope_;}
108  virtual double
109  unprotected_evaluate(IMP::DerivativeAccumulator *accum) const override;
110  virtual IMP::ModelObjectsTemp do_get_inputs() const override;
111  void show(std::ostream &out) const { out << "GEM restraint"; }
112 
113  //! \return Information for writing to RMF files
114  RestraintInfo *get_static_info() const override;
115 
116  //! \return Information for writing to RMF files
117  RestraintInfo *get_dynamic_info() const override;
118 
120 
121  private:
122  double model_cutoff_dist_, density_cutoff_dist_;
123  ParticleIndexes model_ps_;
124  ParticleIndexes density_ps_;
125  Float slope_;
126  bool update_model_,local_;
127  int msize_,dsize_;
128  Float normalization_;
129  Float dd_score_;
130  Float self_mm_score_;
133  ParticleIndexes slope_ps_; //experiment
134  std::string density_fn_;
135 
136  //variables needed to tabulate the exponential
137  Floats exp_grid_;
138  double invdx_;
139  double argmax_;
140 
141  mutable double cross_correlation_;
142 
143  //! Create containers of model and density particles
144  void create_containers();
145 
146  friend class cereal::access;
147 
148  template<class Archive> void serialize(Archive &ar) {
149  ar(cereal::base_class<Restraint>(this), model_cutoff_dist_,
150  density_cutoff_dist_, model_ps_, density_ps_, slope_, update_model_,
151  local_, msize_,dsize_, normalization_, dd_score_, self_mm_score_,
152  slope_ps_, density_fn_, exp_grid_,
153  invdx_, argmax_, cross_correlation_);
154  // recreate md_container and mm_container on input
155  if (std::is_base_of<cereal::detail::InputArchiveBase, Archive>::value) {
156  create_containers();
157  }
158  }
159 
160  IMP_OBJECT_SERIALIZE_DECL(GaussianEMRestraint);
161 };
162 
163 IMPISD_END_NAMESPACE
164 
165 #endif /* IMPISD_GAUSSIAN_EM_RESTRAINT_H */
Decorator to hold Gaussian3D.
A decorator for particles with mass.
virtual RestraintInfo * get_static_info() const
Definition: Restraint.h:178
Store a list of ParticleIndexes.
Return all pairs from a SingletonContainer.
virtual RestraintInfo * get_dynamic_info() const
Definition: Restraint.h:190
#define IMP_OBJECT_METHODS(Name)
Define the basic things needed by any Object.
Definition: object_macros.h:25
Float get_slope()
Get restraint slope.
virtual double unprotected_evaluate(DerivativeAccumulator *da) const
Return the unweighted score for the restraint.
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:86
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.
#define IMP_OBJECT_SERIALIZE_DECL(Name)
Declare methods needed for serialization of Object pointers.
Definition: object_macros.h:95
Report key:value information on restraints.
Definition: RestraintInfo.h:47
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:19
Gaussian shape.
virtual ModelObjectsTemp do_get_inputs() const =0
Class for adding derivatives from restraints to the model.
A restraint is a term in an IMP ScoringFunction.
Definition: Restraint.h:56