IMP logo
IMP Reference Guide  2.12.0
The Integrative Modeling Platform
em_utilities.h
Go to the documentation of this file.
1 /**
2  * \file IMP/isd/em_utilities.h
3  * \brief Common scoring functions
4  *
5  * Copyright 2007-2019 IMP Inventors. All rights reserved.
6  *
7  */
8 
9 #ifndef IMPISD_EM_UTILITIES_H
10 #define IMPISD_EM_UTILITIES_H
11 
12 #include <IMP/algebra/VectorD.h>
13 #include <IMP/algebra/Gaussian3D.h>
15 #include <IMP/core/Gaussian.h>
16 #include <IMP/em/DensityMap.h>
17 #include <IMP/em/DensityHeader.h>
18 #include <boost/random/uniform_real.hpp>
19 #include <boost/random/variate_generator.hpp>
20 
21 IMPISD_BEGIN_NAMESPACE
22 
23 #if !defined(SWIG) && !defined(IMP_DOXYGEN)
24 inline Float score_gaussian_overlap(Model *m,
25  ParticleIndexPair pp,
26  Eigen::Vector3d * deriv){
27  double determinant;
28  bool invertible;
29  Eigen::Matrix3d inverse = Eigen::Matrix3d::Zero();
30  Float mass12 = atom::Mass(m,pp[0]).get_mass() *
31  atom::Mass(m,pp[1]).get_mass();
32  core::Gaussian g1(m,pp[0]);
33  core::Gaussian g2(m,pp[1]);
34  Eigen::Matrix3d covar = g1.get_global_covariance() +
35  g2.get_global_covariance();
36  Eigen::Vector3d v = Eigen::Vector3d(g2.get_coordinates().get_data())
37  - Eigen::Vector3d(g1.get_coordinates().get_data());
38  covar.computeInverseAndDetWithCheck(inverse,determinant,invertible);
39  Eigen::Vector3d tmp = inverse*v;
40  // 0.06349... = 1. / sqrt(2.0 * pi) ** 3
41  Float score = mass12 * 0.06349363593424097 / (std::sqrt(determinant)) *
42  std::exp(-0.5*v.transpose()*tmp);
43  *deriv = -score*tmp;
44  return score;
45 }
46 #endif
47 
48 inline FloatsList sample_points_from_density(const em::DensityMap * dmap_orig,
49  int npoints,
50  Float threshold=0.0){
51  // get sample region
52  em::DensityMap * dmap = em::get_threshold_map(dmap_orig,threshold);
53  dmap->calcRMS();
54  const em::DensityHeader * dhead = dmap->get_header();
55  Float dmax=dhead->dmax;
56  algebra::BoundingBox3D bbox=em::get_bounding_box(dmap,0.00001);
57 
58  // setup random number generator
59  FloatsList ret;
60  boost::mt19937 generator(std::time(0));
61  boost::uniform_real<> uni_dist(0,1);
62  boost::variate_generator<boost::mt19937&, boost::uniform_real<> > uni(generator, uni_dist);
63  for (int i=0;i<npoints;i++){
65  Float den=(dmap->get_value(vs))/dmax;
66  Float t=uni();
67  if (t<den) {
68  Floats r;
69  r.push_back(vs[0]);
70  r.push_back(vs[1]);
71  r.push_back(vs[2]);
72  ret.push_back(r);
73  }
74  }
75  return ret;
76 }
77 
78 IMPISD_END_NAMESPACE
79 
80 #endif /* IMPISD_EM_UTILITIES_H */
Decorator to hold Gaussian3D.
BoundingBoxD< 3 > BoundingBox3D
Typedef for Python.
Definition: BoundingBoxD.h:176
IMP::Vector< Float > Floats
Standard way to pass a bunch of Float values.
Definition: types.h:47
Metadata for a density file.
DensityMap * get_threshold_map(const DensityMap *orig_map, float threshold)
Return a map with 0 for all voxels below the threshold.
Class for handling density maps.
IMP::Vector< Floats > FloatsList
Standard way to pass a bunch of Floats values.
Definition: types.h:54
Functions to generate vectors.
VectorD< D > get_random_vector_in(const BoundingBoxD< D > &bb)
Generate a random vector in a box with uniform density.
Simple D vector class.
VectorD< 3 > Vector3D
Definition: VectorD.h:421
double Float
Basic floating-point value (could be float, double...)
Definition: types.h:20
Gaussian shape.