IMP logo
IMP Reference Guide  2.6.1
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-2016 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>
14 #include <IMP/core/Gaussian.h>
15 #include <IMP/em/DensityMap.h>
16 #include <IMP/em/DensityHeader.h>
17 #include <boost/random/uniform_real.hpp>
18 #include <boost/random/variate_generator.hpp>
19 
20 IMPISD_BEGIN_NAMESPACE
21 
22 #if !defined(SWIG) && !defined(IMP_DOXYGEN)
23 inline Float score_gaussian_overlap(Model *m,
24  ParticleIndexPair pp,
25  IMP_Eigen::Vector3d * deriv){
26  double determinant;
27  bool invertible;
28  IMP_Eigen::Matrix3d inverse = IMP_Eigen::Matrix3d::Zero();
29  Float mass12 = atom::Mass(m,pp[0]).get_mass() *
30  atom::Mass(m,pp[1]).get_mass();
31  core::Gaussian g1(m,pp[0]);
32  core::Gaussian g2(m,pp[1]);
33  IMP_Eigen::Matrix3d covar = g1.get_global_covariance() +
34  g2.get_global_covariance();
35  IMP_Eigen::Vector3d v = IMP_Eigen::Vector3d(g2.get_coordinates().get_data())
36  - IMP_Eigen::Vector3d(g1.get_coordinates().get_data());
37  covar.computeInverseAndDetWithCheck(inverse,determinant,invertible);
38  IMP_Eigen::Vector3d tmp = inverse*v;
39  Float score = mass12 * 0.06349363593424097 / (std::sqrt(determinant)) *
40  std::exp(-0.5*v.transpose()*tmp);
41  *deriv = -score*tmp;
42  return score;
43 }
44 #endif
45 
46 inline FloatsList sample_points_from_density(const em::DensityMap * dmap_orig,
47  int npoints,
48  Float threshold=0.0){
49  // get sample region
50  em::DensityMap * dmap = em::get_threshold_map(dmap_orig,threshold);
51  dmap->calcRMS();
52  const em::DensityHeader * dhead = dmap->get_header();
53  Float dmax=dhead->dmax;
54  algebra::BoundingBox3D bbox=em::get_bounding_box(dmap,0.00001);
55 
56  // setup random number generator
57  FloatsList ret;
58  boost::mt19937 generator(std::time(0));
59  boost::uniform_real<> uni_dist(0,1);
60  boost::variate_generator<boost::mt19937&, boost::uniform_real<> > uni(generator, uni_dist);
61  for (int i=0;i<npoints;i++){
63  Float den=(dmap->get_value(vs))/dmax;
64  Float t=uni();
65  if (t<den) {
66  Floats r;
67  r.push_back(vs[0]);
68  r.push_back(vs[1]);
69  r.push_back(vs[2]);
70  ret.push_back(r);
71  }
72  }
73  return ret;
74 }
75 
76 IMPISD_END_NAMESPACE
77 
78 #endif /* IMPISD_EM_UTILITIES_H */
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
VectorD< D > get_random_vector_in(const BoundingBoxD< D > &bb)
Generate a random vector in a box with uniform density.
Simple D vector class.
Decorator to hold Gaussian3D.
VectorD< 3 > Vector3D
Definition: VectorD.h:395
double Float
Basic floating-point value (could be float, double...)
Definition: types.h:20
Gaussian shape.