IMP logo
IMP Reference Guide  2.6.0
The Integrative Modeling Platform
GaussianEMRestraint.h
Go to the documentation of this file.
1 /**
2  * \file IMP/isd/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 IMPISD_GAUSSIAN_EM_RESTRAINT_H
10 #define IMPISD_GAUSSIAN_EM_RESTRAINT_H
11 
12 #include "isd_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 <IMP/algebra/eigen3/Eigen/Dense>
24 #include <boost/unordered_map.hpp>
25 
26 IMPISD_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  IMP_Eigen::Vector3d sum;
38  IMP_Eigen::Vector3d correction;
39 KahanVectorAccumulation():
40  sum(IMP_Eigen::Vector3d(0,0,0)),
41  correction(IMP_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, IMP_Eigen::Vector3d value){
54  KahanVectorAccumulation result;
55  IMP_Eigen::Vector3d y = value - accumulation.correction;
56  IMP_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 IMPISDEXPORT 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  //! Pre-calculate the density-density and model-model scores
108  /** This is automatically called by the constructor.
109  You only need to call it manually if you change Gaussian variances
110  */
111  void compute_initial_scores();
112 
113  //! Set restraint slope
114  void set_slope(Float s){slope_=s;}
115 
116  //! Get restraint slope
117  Float get_slope(){return slope_;}
118  virtual double
119  unprotected_evaluate(IMP::DerivativeAccumulator *accum)
120  const IMP_OVERRIDE;
122  void show(std::ostream &out) const { out << "GEM restraint"; }
124 
125  private:
126  ParticleIndexes model_ps_;
127  ParticleIndexes density_ps_;
128  ParticleIndex global_sigma_;
129  Float slope_;
130  bool update_model_;
131  int msize_,dsize_;
132  Float normalization_;
133  Float dd_score_;
134  Float self_mm_score_;
137  ParticleIndexes slope_ps_; //experiment
138 
139  //variables needed to tabulate the exponential
140  Floats exp_grid_;
141  double invdx_;
142  double argmax_;
143 
144 };
145 
146 IMPISD_END_NAMESPACE
147 
148 #endif /* IMPISD_GAUSSIAN_EM_RESTRAINT_H */
Declare an efficient stl-compatible map.
A decorator for particles with mass.
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_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.
Creates a restraint between two Gaussian Mixture Models, "model" and "density".
Return all pairs from a SingletonContainer.
Decorator to hold Gaussian3D.
std::ostream & show(Hierarchy h, std::ostream &out=std::cout)
Print the hierarchy using a given decorator to display each node.
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:52