IMP logo
IMP Reference Guide  develop.cb6747d2d1,2024/03/28
The Integrative Modeling Platform
bayesianem/GaussianEMRestraint.h
Go to the documentation of this file.
1 /**
2  * \file IMP/bayesianem/GaussianEMRestraint.h
3  * \brief Restraint two sets of gaussians (model and gmm derived from EM map)
4  *
5  * Copyright 2007-2016 IMP Inventors. All rights reserved.
6  *
7  */
8 
9 #ifndef IMPBAYESIANEM_GAUSSIAN_EM_RESTRAINT_H
10 #define IMPBAYESIANEM_GAUSSIAN_EM_RESTRAINT_H
11 
12 #include "bayesianem_config.h"
13 #include <IMP/PairContainer.h>
15 #include <IMP/container_macros.h>
16 #include <IMP/core/XYZ.h>
17 #include <IMP/core/Gaussian.h>
18 #include <IMP/algebra/Gaussian3D.h>
21 #include <IMP/atom/Mass.h>
22 #include <math.h>
23 #include <Eigen/Dense>
24 #include <boost/unordered_map.hpp>
25 
26 IMPBAYESIANEM_BEGIN_NAMESPACE
27 
28 #if !defined(SWIG) && !defined(IMP_DOXYGEN)
29 struct KahanAccumulation{
30 double sum;
31 double correction;
32 KahanAccumulation():
33  sum(0.0),correction(0.0)
34  {}
35 };
36 struct KahanVectorAccumulation{
37  Eigen::Vector3d sum;
38  Eigen::Vector3d correction;
39 KahanVectorAccumulation():
40  sum(Eigen::Vector3d(0,0,0)),
41  correction(Eigen::Vector3d(0,0,0))
42  {}
43 };
44 inline KahanAccumulation KahanSum(KahanAccumulation accumulation, double value){
45  KahanAccumulation result;
46  double y = value - accumulation.correction;
47  double t = accumulation.sum + y;
48  result.correction = (t - accumulation.sum) - y;
49  result.sum = t;
50  return result;
51 }
52 inline KahanVectorAccumulation
53 KahanVectorSum(KahanVectorAccumulation accumulation, Eigen::Vector3d value){
54  KahanVectorAccumulation result;
55  Eigen::Vector3d y = value - accumulation.correction;
56  Eigen::Vector3d t = accumulation.sum + y;
57  result.correction = (t - accumulation.sum) - y;
58  result.sum = t;
59  return result;
60 }
61 #endif
62 
63 
64 //! Creates a restraint between two Gaussian Mixture Models, "model" and "density"
65 //
66 /** This restrains two sets of GMMs with a density overlap function.
67  The function is correlation of the two GMMs \f$f_M\f$ and \f$f_D\f$:
68  \f[
69  \frac{2\int{f_M(x)f_D(x)dx}}{\int{f_M^2(x)+f_D^2(x)dx}}
70  \f]
71  Where the integral is the "overlap function" given by:
72  \f[
73  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 ]
74  \f]
75  \note Source: Greenberg, Pellarin, Sali. In preparation.
76  */
77 class IMPBAYESIANEMEXPORT GaussianEMRestraint : public Restraint
78 {
79 
80  public:
81  //! Setup the GaussianEMRestraint
82  /**
83  \param[in] mdl the Model object to operate on
84  \param[in] model_ps particles for the model GMM
85  \param[in] density_ps particles for the density GMM
86  \param[in] global_sigma Particle to modulate the uncertainty
87  \param[in] model_cutoff_dist Cutoff for the model-model interactions
88  \param[in] density_cutoff_dist Cutoff for model-density interactions
89  \param[in] slope Gentle term to move all particles to the density
90  \param[in] update_model (DEPRECATED) update model each cycle
91  \param[in] backbone_slope Limit the slope only to backbone particles
92  \param[in] name Name of this restraint
93  \note the model and density particles must be set up as Gaussian
94  */
96  ParticleIndexes model_ps, ParticleIndexes density_ps,
97  ParticleIndex global_sigma,
98  Float model_cutoff_dist,Float density_cutoff_dist,Float slope,
99  bool update_model=true, bool backbone_slope=false,
100  std::string name="GaussianEMRestraint%1%");
101 
102  //! Returns exp(score)
103  double get_probability() const {
104  return exp(-unprotected_evaluate(NULL));
105  }
106 
107  //! Set the filename corresponding to the density GMM particles
108  /** If the density GMM particles were read from a file, this method
109  can be used to tell the restraint so that it can track this
110  information back to the original EM density file, which is useful
111  for deposition. */
112  void set_density_filename(std::string density_fn) {
113  density_fn_ = get_absolute_path(density_fn);
114  }
115 
116  //! Pre-calculate the density-density and model-model scores
117  /** This is automatically called by the constructor.
118  You only need to call it manually if you change Gaussian variances
119  */
120  void compute_initial_scores();
121 
122  ParticleIndexes const get_indexes(){
123  return pis_debug_;
124  };
125  Floats const get_log2(){
126  return log2_debug_;
127  };
128 
129  void debug();
130  //! Set restraint slope
131  void set_slope(Float s){slope_=s;}
132 
133  //! Get restraint slope
134  Float get_slope(){return slope_;}
135 
136  virtual double
137  unprotected_evaluate(IMP::DerivativeAccumulator *accum) const override;
138  virtual IMP::ModelObjectsTemp do_get_inputs() const override;
139  void show(std::ostream &out) const { out << "GEM restraint"; }
140 
141  //! \return Information for writing to RMF files
142  RestraintInfo *get_static_info() const override;
143 
144  //! \return Information for writing to RMF files
145  RestraintInfo *get_dynamic_info() const override;
146 
148 
149  private:
150  double model_cutoff_dist_, density_cutoff_dist_;
151  ParticleIndexes model_ps_;
152  ParticleIndexes density_ps_;
153  ParticleIndex global_sigma_;
154  ParticleIndexes pis_debug_;
155  Floats log2_debug_;
156  Float slope_;
157  bool update_model_;
158  int msize_,dsize_;
159  Float normalization_;
160  Float dd_score_;
161  Float self_mm_score_;
164  ParticleIndexes slope_ps_; //experiment
165  std::map<ParticleIndex,Float> map_score_dd_;
166  Float cached_score_term_;
167  std::string density_fn_;
168 
169  //variables needed to tabulate the exponential
170  Floats exp_grid_;
171  double invdx_;
172  double argmax_;
173 
174 
175 };
176 
177 IMPBAYESIANEM_END_NAMESPACE
178 
179 #endif /* IMPBAYESIANEM_GAUSSIAN_EM_RESTRAINT_H */
Helper functions to check for NaN or infinity.
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
Creates a restraint between two Gaussian Mixture Models, "model" and "density".
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.
Macros to define containers of objects.
Simple XYZ decorator.
Class for storing model, its restraints, constraints, and particles.
Definition: Model.h:86
A container for Pairs.
Return all pairs from a SingletonContainer.
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.
void set_slope(Float s)
Set restraint slope.
double get_probability() const
Returns exp(score)
ParticleIndexes get_indexes(const ParticlesTemp &ps)
Get the indexes from a list of particles.
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