IMP logo
IMP Reference Guide  2.18.0
The Integrative Modeling Platform
GaussianAnchorEMRestraint.h
Go to the documentation of this file.
1 /**
2  * \file IMP/isd/GaussianAnchorEMRestraint.h
3  * \brief Restraint between two sets of anchor points "model" and "EM density"
4  *
5  * Copyright 2007-2022 IMP Inventors. All rights reserved.
6  *
7  */
8 
9 #ifndef IMPISD_GAUSSIAN_ANCHOR_EMRESTRAINT_H
10 #define IMPISD_GAUSSIAN_ANCHOR_EMRESTRAINT_H
11 
12 #include <IMP/isd/isd_config.h>
13 #include <IMP/isd/ISDRestraint.h>
14 #include <IMP/PairContainer.h>
15 #include <IMP/isd/Scale.h>
16 #include <IMP/core/XYZ.h>
18 #include <map>
19 #include <math.h>
20 
21 IMPISD_BEGIN_NAMESPACE
22 
23 //! Restraint between two sets of anchor points "model" and "EM density"
24 /** \note This class is obsolete - use GaussianEMRestraint instead. It exists
25  solely to support older applications of IMP that used this method.
26  */
27 class IMPISDEXPORT GaussianAnchorEMRestraint : public ISDRestraint
28 {
29 public:
30  //! Constructor
31  /**
32  \param[in] model_ps particles for the model GMM
33  \param[in] model_sigs sigmas for the model GMM
34  (currently spherically symmetric)
35  \param[in] model_weights weights for the model GMM
36  \param[in] density_ps particles for the density GMM
37  \param[in] density_sigs sigmas for the density GMM
38  (currently spherically symmetric)
39  \param[in] density_weights weights for the density GMM
40  \param[in] global_sigma Particle to modulate the uncertainty
41  \param[in] rigid Set to true if the model is rigid (faster)
42  */
44  Floats model_weights, ParticlesTemp density_ps,
45  Floats density_sigs, Floats density_weights,
46  Particle *global_sigma,Float cutoff_dist,
47  bool rigid, bool tabexp,
48  std::string name="GaussianAnchorEMRestraint%1%"):
49  ISDRestraint(model_ps[0]->get_model(),name),
50  model_ps_(model_ps),
51  model_sigs_(model_sigs),
52  model_weights_(model_weights),
53  density_ps_(density_ps),
54  density_sigs_(density_sigs),
55  density_weights_(density_weights),
56  global_sigma_(global_sigma),
57  cutoff_dist_(cutoff_dist),
58  rigid_(rigid),
59  tabexp_(tabexp){
60  msize_=model_ps.size();
61  dsize_=density_ps.size();
62  IMP_USAGE_CHECK(model_sigs.size()==msize_ &&
63  model_weights.size()==msize_,
64  "All model inputs must be same size");
65  IMP_USAGE_CHECK(density_ps.size()==dsize_ &&
66  density_sigs.size()==dsize_ &&
67  density_weights.size()==dsize_,
68  "All input vectors must be same size");
69 
70  if (tabexp_){
71  unsigned exparg_grid_size=1000001;
72  argmax_=10.0;
73  invdx_=double(exparg_grid_size)/argmax_;
74  for(unsigned k=0;k<exparg_grid_size;++k){
75  double argvalue=double(k)/invdx_;
76  exp_grid_.push_back(std::exp(-argvalue));
77  }
78  }
79 
80 
81  // calculating prefactors and initial model-model score
82  init_mm_score_=0.0;
83  for (size_t m1=0;m1<msize_;m1++){
84  init_mm_score_+=calc_prefactor(model_sigs_[m1],model_weights_[m1],
85  model_sigs[m1],model_weights_[m1]);
86  for (size_t m2=m1;m2<msize_;m2++){
87  Float dist=core::get_distance(core::XYZ(model_ps_[m1]),
88  core::XYZ(model_ps_[m2]));
89  ParticlePair pp1(model_ps_[m1],model_ps_[m2]);
90  ParticlePair pp2(model_ps_[m2],model_ps_[m1]);
91  Float pref=calc_prefactor(model_sigs_[m1],model_weights_[m1],
92  model_sigs[m2],model_weights_[m2]);
93  Float prod=calc_prod(model_sigs_[m1],model_sigs_[m2]);
94  mm_prefactors_[pp1]=pref;
95  mm_prefactors_[pp2]=pref;
96  mm_prods_[pp1]=prod;
97  mm_prods_[pp2]=prod;
98  if (m2>m1){
99  init_mm_score_+=2.0*pref*calc_score(dist,prod);
100  }
101  }
102 
103  //model-density prefactors
104  for (size_t d1=0;d1<dsize_;d1++){
105  ParticlePair ppd1(model_ps_[m1],density_ps_[d1]);
106  md_prefactors_[ppd1]=
107  calc_prefactor(model_sigs_[m1],model_weights_[m1],
108  density_sigs[d1],density_weights_[d1]);
109  md_prods_[ppd1]=calc_prod(model_sigs_[m1],density_sigs[d1]);
110  }
111  }
112 
113  // calculating density-density score (doesn't change)
114  for (size_t d1=0;d1<dsize_;d1++){
115  core::XYZ d_d1(density_ps_[d1]);
116  dd_score_+=calc_prefactor(density_sigs_[d1],density_weights_[d1],
117  density_sigs_[d1],density_weights_[d1]);
118  for (size_t d2=d1+1;d2<dsize_;d2++){
119  Float dist=core::get_distance(core::XYZ(density_ps_[d1]),
120  core::XYZ(density_ps_[d2]));
121  Float pref=calc_prefactor(density_sigs_[d1],density_weights_[d1],
122  density_sigs_[d2],density_weights_[d2]);
123  Float prod=calc_prod(density_sigs_[d1],density_sigs_[d2]);
124  dd_score_+=2.0*pref*calc_score(dist,prod);
125  }
126  }
127  //KLUGE!
128  if (rigid_){
129  init_mm_score_=0.0;
130  dd_score_=0.0;
131  }
132 
133  //Set up md container
134  md_container_ = new container::CloseBipartitePairContainer(
135  model_ps,density_ps,cutoff_dist_);
136  }
137  ParticlesTemp get_density_particles() const {
138  ParticlesTemp ret = density_ps_;
139  return ret;
140  }
141  double get_probability() const override {
142  return exp(-unprotected_evaluate(NULL));
143  }
144 
145  virtual double
146  unprotected_evaluate(IMP::DerivativeAccumulator *accum) const override;
147  virtual IMP::ModelObjectsTemp do_get_inputs() const override;
149  protected:
150  Particles model_ps_;
151  Floats model_sigs_,model_weights_;
152  Particles density_ps_;
153  Floats density_sigs_,density_weights_;
154  Particle * global_sigma_;
155  Float cutoff_dist_;
156  bool rigid_;
157  bool tabexp_;
158  Float dd_score_,init_mm_score_;
159  std::map<ParticlePair,Float> mm_prefactors_,
160  md_prefactors_,mm_prods_,md_prods_;
161  size_t msize_,dsize_;
162  //variables needed to tabulate the exponential
163  Floats exp_grid_;
164  double invdx_;
165  double argmax_;
167  //actual scores defined here (prefactor*score)
168  inline Float calc_prefactor(Float s1, Float w1, Float s2, Float w2) const{
169  return w1*w2*sqrt(8*algebra::PI)*s1*s1*s2*s2/
170  pow(s1*s1+s2*s2,1.5);
171  }
172  inline Float calc_prod(Float s1, Float s2) const{
173  return 1.0/(2.0*(s1*s1+s2*s2));
174  }
175  inline Float calc_score (Float dist, Float prod) const{
176  double argvalue=dist*dist*prod;
177  double expn;
178  if (tabexp_){
179  double minarg=std::min(argvalue,argmax_);
180  unsigned k = static_cast<unsigned>( std::floor(minarg*invdx_) );
181  expn=exp_grid_[k];
182  }
183  else{
184  expn=std::exp(-argvalue);
185  }
186  return expn;
187  }
188 };
189 
190 IMPISD_END_NAMESPACE
191 
192 #endif /* IMPISD_GAUSSIAN_ANCHOR_EMRESTRAINT_H */
Helper functions to check for NaN or infinity.
A base class for ISD Restraints.
static const double PI
the constant pi
A class to store an fixed array of same-typed values.
Definition: Array.h:33
#define IMP_OBJECT_METHODS(Name)
Define the basic things needed by any Object.
Definition: object_macros.h:25
GaussianAnchorEMRestraint(ParticlesTemp model_ps, Floats model_sigs, Floats model_weights, ParticlesTemp density_ps, Floats density_sigs, Floats density_weights, Particle *global_sigma, Float cutoff_dist, bool rigid, bool tabexp, std::string name="GaussianAnchorEMRestraint%1%")
Constructor.
A decorator for scale parameters particles.
Return all spatially-proximals pairs of particles (a,b) from the two SingletonContainers A and B...
Simple XYZ decorator.
virtual ModelObjectsTemp do_get_inputs() const override
A container for Pairs.
Return all pairs from a SingletonContainer.
Restraint between two sets of anchor points "model" and "EM density".
A decorator for a particle with x,y,z coordinates.
Definition: XYZ.h:30
virtual double unprotected_evaluate(DerivativeAccumulator *accum) const override
Return the unweighted score for the restraint.
double Float
Basic floating-point value (could be float, double...)
Definition: types.h:19
Class to handle individual particles of a Model object.
Definition: Particle.h:41
#define IMP_USAGE_CHECK(expr, message)
A runtime test for incorrect usage of a class or method.
Definition: check_macros.h:168
double get_distance(const Surface &s, const XYZR &d)
Get distance from sphere to surface.
Definition: Surface.h:153
Class for adding derivatives from restraints to the model.