IMP logo
IMP Reference Guide  2.14.0
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)
138  const IMP_OVERRIDE;
140  void show(std::ostream &out) const { out << "GEM restraint"; }
141 
142  //! \return Information for writing to RMF files
144 
145  //! \return Information for writing to RMF files
146  RestraintInfo *get_dynamic_info() const IMP_OVERRIDE;
147 
149 
150  private:
151  double model_cutoff_dist_, density_cutoff_dist_;
152  ParticleIndexes model_ps_;
153  ParticleIndexes density_ps_;
154  ParticleIndex global_sigma_;
155  ParticleIndexes pis_debug_;
156  Floats log2_debug_;
157  Float slope_;
158  bool update_model_;
159  int msize_,dsize_;
160  Float normalization_;
161  Float dd_score_;
162  Float self_mm_score_;
163  PointerMember<container::CloseBipartitePairContainer> md_container_;
164  PointerMember<container::ClosePairContainer> mm_container_;
165  ParticleIndexes slope_ps_; //experiment
166  std::map<ParticleIndex,Float> map_score_dd_;
167  Float cached_score_term_;
168  std::string density_fn_;
169 
170  //variables needed to tabulate the exponential
171  Floats exp_grid_;
172  double invdx_;
173  double argmax_;
174 
175 
176 };
177 
178 IMPBAYESIANEM_END_NAMESPACE
179 
180 #endif /* IMPBAYESIANEM_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
Creates a restraint between two Gaussian Mixture Models, "model" and "density".
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:72
A container for Pairs.
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.
void set_slope(Float s)
Set restraint slope.
double get_probability() const
Returns exp(score)
ParticleIndexes get_indexes(const ParticlesTemp &ps)
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