IMP logo
IMP Reference Guide  develop.63b38c487d,2024/12/22
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-2022 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 <IMP/random.h>
19 #include <boost/random/uniform_real_distribution.hpp>
20 #include <boost/random/variate_generator.hpp>
21 
22 IMPISD_BEGIN_NAMESPACE
23 
24 #if !defined(SWIG) && !defined(IMP_DOXYGEN)
25 inline Float score_gaussian_overlap(Model *m,
26  ParticleIndexPair pp,
27  Eigen::Vector3d * deriv){
28  double determinant;
29  bool invertible;
30  Eigen::Matrix3d inverse = Eigen::Matrix3d::Zero();
31  Float mass12 = atom::Mass(m, std::get<0>(pp)).get_mass() *
32  atom::Mass(m, std::get<1>(pp)).get_mass();
33  core::Gaussian g1(m, std::get<0>(pp));
34  core::Gaussian g2(m, std::get<1>(pp));
35  Eigen::Matrix3d covar = g1.get_global_covariance() +
36  g2.get_global_covariance();
37  Eigen::Vector3d v = Eigen::Vector3d(g2.get_coordinates().get_data())
38  - Eigen::Vector3d(g1.get_coordinates().get_data());
39  covar.computeInverseAndDetWithCheck(inverse,determinant,invertible);
40  Eigen::Vector3d tmp = inverse*v;
41  // 0.06349... = 1. / sqrt(2.0 * pi) ** 3
42  Float score = mass12 * 0.06349363593424097 / (std::sqrt(determinant)) *
43  std::exp(-0.5*v.transpose()*tmp);
44  *deriv = -score*tmp;
45  return score;
46 }
47 #endif
48 
49 inline FloatsList sample_points_from_density(const em::DensityMap * dmap_orig,
50  int npoints,
51  Float threshold=0.0){
52  // get sample region
53  em::DensityMap * dmap = em::get_threshold_map(dmap_orig,threshold);
54  dmap->calcRMS();
55  const em::DensityHeader * dhead = dmap->get_header();
56  Float dmax=dhead->dmax;
57  algebra::BoundingBox3D bbox=em::get_bounding_box(dmap,0.00001);
58 
59  // setup random number generator
60  FloatsList ret;
61  boost::random::uniform_real_distribution<> uni_dist(0,1);
62  boost::variate_generator<
63  IMP::RandomNumberGenerator&,
64  boost::random::uniform_real_distribution<> >
65  uni(IMP::random_number_generator, uni_dist);
66  for (int i=0;i<npoints;i++){
68  Float den=(dmap->get_value(vs))/dmax;
69  Float t=uni();
70  if (t<den) {
71  Floats r;
72  r.push_back(vs[0]);
73  r.push_back(vs[1]);
74  r.push_back(vs[2]);
75  ret.push_back(r);
76  }
77  }
78  return ret;
79 }
80 
81 IMPISD_END_NAMESPACE
82 
83 #endif /* IMPISD_EM_UTILITIES_H */
Decorator to hold Gaussian3D.
BoundingBoxD< 3 > BoundingBox3D
Typedef for Python.
Definition: BoundingBoxD.h:183
IMP::Vector< Float > Floats
Standard way to pass a bunch of Float values.
Definition: types.h:46
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:53
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:408
double Float
Basic floating-point value (could be float, double...)
Definition: types.h:19
Random number generators used by IMP.
Gaussian shape.
RandomNumberGenerator random_number_generator
A shared non-GPU random number generator.