IMP logo
IMP Reference Guide  2.13.0
The Integrative Modeling Platform
bayesianem/utilities.h
Go to the documentation of this file.
1 /**
2  * \file IMP/bayesianem/utilities.h
3  * \brief Useful utilities
4  *
5  * Copyright 2007-2016 IMP Inventors. All rights reserved.
6  */
7 
8 #ifndef IMPBAYESIANEM_UTILITIES_H
9 #define IMPBAYESIANEM_UTILITIES_H
10 
11 #include <IMP/bayesianem/bayesianem_config.h>
12 //#include <IMP/algebra/algebra_config.h>
13 #include "IMP/algebra/Gaussian3D.h"
15 #include <Eigen/Geometry>
16 #include <Eigen/LU>
17 #include <Eigen/Eigenvalues>
18 #include <IMP/em.h>
19 #include <IMP/isd.h>
20 
21 #include <limits>
22 #include <vector>
23 
24 IMPBAYESIANEM_BEGIN_NAMESPACE
25 
26 typedef IMP::algebra::DenseGrid3D<double> DensityGrid;
27 
28 double get_rmsd_excluding_higher_than( const IMP::core::XYZs &m1,
29  const IMP::core::XYZs &m2,
30  const double t) {
31  double rmsd = 0.0;
32  int N(0);
33  for (int i=0; i<m1.size(); ++i) {
34  double tmprmsd = get_squared_distance(m1[i].get_coordinates(), m2[i].get_coordinates());
35  if(tmprmsd<t*t){
36  rmsd+=tmprmsd;
37  ++N;
38  }
39  }
40  return std::sqrt(rmsd / N);
41  //return Floats({std::sqrt(rmsd / N),((double)N/m1.size())});
42 }
43 
44 double get_percentage_closer_than(const IMP::core::XYZs &m1,
45  const IMP::core::XYZs &m2,
46  const double t) {
47  int N(0);
48  for (int i=0; i<m1.size(); ++i) {
49  double tmprmsd = get_squared_distance(m1[i].get_coordinates(), m2[i].get_coordinates());
50  if(tmprmsd<t*t){
51  ++N;
52  }
53  }
54  return ((double) N)/m1.size();
55  //return Floats({std::sqrt(rmsd / N),((double)N/m1.size())});
56 }
57 
58 double get_rmsd_of_best_population(const IMP::core::XYZs &m1,
59  const IMP::core::XYZs &m2,
60  const double percentage){
61  std::vector<double> sq_distances(m1.size());
62  for (int i=0; i<m1.size(); ++i) {
63  sq_distances[i] = get_squared_distance(m1[i].get_coordinates(), m2[i].get_coordinates());
64  }
65  std::sort(sq_distances.begin(), sq_distances.end());
66  double sd=0.0;
67  int N=0;
68  for (int i=0; ((double) i) < (percentage*m1.size()); ++i) {
69  sd+=sq_distances[i];
70  ++N;
71  }
72  return std::sqrt(sd/N);
73 }
74 
75 double get_rmsd_of_best_population(const IMP::atom::Selection &s1,
76  const IMP::atom::Selection &s2,
77  const double percentage){
78 
81 
82  IMP::core::XYZs ds1;
83  IMP::core::XYZs ds2;
84 
85  for (int i=0; i<ps1.size(); ++i){ds1.push_back(IMP::core::XYZ(ps1[i]));}
86  for (int i=0; i<ps2.size(); ++i){ds2.push_back(IMP::core::XYZ(ps2[i]));}
87 
88  return get_rmsd_of_best_population(ds1,ds2,percentage);
89 }
90 
92  const IMP::atom::Selection &s2){
93 
96 
97  IMP::core::XYZs ds1;
98  IMP::core::XYZs ds2;
99 
100  for (int i=0; i<ps1.size(); ++i){ds1.push_back(IMP::core::XYZ(ps1[i]));}
101  for (int i=0; i<ps2.size(); ++i){ds2.push_back(IMP::core::XYZ(ps2[i]));}
102 
104 }
105 
106 
107 
108 double get_rmsd_of_best_population(const IMP::algebra::Vector3Ds &m1,
109  const IMP::algebra::Vector3Ds &m2,
110  const double percentage){
111  std::vector<double> sq_distances(m1.size());
112  for (int i=0; i<m1.size(); ++i) {
113  sq_distances[i] = get_squared_distance(m1[i], m2[i]);
114  }
115  std::sort(sq_distances.begin(), sq_distances.end());
116  double sd=0.0;
117  int N=0;
118  for (int i=0; ((double) i) < (percentage*m1.size()); ++i) {
119  sd+=sq_distances[i];
120  ++N;
121  }
122  return std::sqrt(sd/N);
123 }
124 
125 double gem_score_cc(Particles model_ps, Particles density_ps){
126  Eigen::Vector3d deriv;
127  double mm_score(0.0);
128  double md_score(0.0);
129  double dd_score(0.0);
130 
131  int nm = model_ps.size();
132  int nd = density_ps.size();
133  IMP::Model *mdl = model_ps[0]->get_model();
134 
135  for(int nd1=0 ; nd1<nd ; ++nd1){
136  for(int nd2=0 ; nd2<nd ; ++nd2){
137  dd_score += IMP::isd::score_gaussian_overlap(mdl, ParticleIndexPair(density_ps[nd1]->get_index(), density_ps[nd2]->get_index()), &deriv);
138  }
139  }
140  for(int nm1=0 ; nm1<nm ; ++nm1){
141  for(int nm2=0 ; nm2<nm ; ++nm2){
142  mm_score += IMP::isd::score_gaussian_overlap(mdl, ParticleIndexPair(model_ps[nm1]->get_index(), model_ps[nm2]->get_index()), &deriv);
143  }
144  for(int nd1=0 ; nd1<nd ; ++nd1){
145  md_score += IMP::isd::score_gaussian_overlap(mdl, ParticleIndexPair(model_ps[nm1]->get_index(), density_ps[nd1]->get_index()), &deriv);
146  }
147  }
148  double cc = md_score/std::sqrt(mm_score*dd_score);
149  return cc;
150 }
151 
152 
153 namespace {
154 double get_gaussian_eval_prefactor(double determinant) {
155  return 1.0 / pow(2 * algebra::PI, 2.0 / 3.0) / std::sqrt(determinant);
156 }
157 Eigen::Vector3d get_vec_from_center(IMP::algebra::GridIndex3D i,
158  DensityGrid const &g,
159  Eigen::Vector3d const &center) {
160  IMP::algebra::Vector3D aloc = g.get_center(i);
161  Eigen::Vector3d loc(aloc[0], aloc[1], aloc[2]);
162  Eigen::Vector3d r(loc - center);
163  return r;
164 }
165 void update_value(DensityGrid *g,
166  DensityGrid::Index i, Eigen::Vector3d r,
167  Eigen::Matrix3d inverse, double pre, Float weight) {
168  double d(r.transpose() * (inverse * r));
169  double score(pre * weight * std::exp(-0.5 * (d)));
170  (*g)[i] += score;
171 }
172 }
173 
174 DensityGrid get_grid(IMP::em::DensityMap *in) {
176  IMP_CHECK_OBJECT(in);
177  DensityGrid ret(in->get_header()->get_spacing(), get_bounding_box(in));
178  IMP_USAGE_CHECK(ret.get_number_of_voxels(0) ==
179  static_cast<unsigned int>(in->get_header()->get_nx()),
180  "X voxels don't match");
181  IMP_USAGE_CHECK(ret.get_number_of_voxels(1) ==
182  static_cast<unsigned int>(in->get_header()->get_ny()),
183  "Y voxels don't match");
184  IMP_USAGE_CHECK(ret.get_number_of_voxels(2) ==
185  static_cast<unsigned int>(in->get_header()->get_nz()),
186  "Z voxels don't match");
187  for (unsigned int i = 0; i < ret.get_number_of_voxels(0); ++i) {
188  for (unsigned int j = 0; j < ret.get_number_of_voxels(1); ++j) {
189  for (unsigned int k = 0; k < ret.get_number_of_voxels(2); ++k) {
190  DensityGrid::ExtendedIndex ei(i, j, k);
191  DensityGrid::Index gi = ret.get_index(ei);
192  long vi = in->get_voxel_by_location(i, j, k);
193  ret[gi] = in->get_value(vi);
194  }
195  }
196  }
197  return ret;
198 }
199 
200 IMP::em::DensityMap *get_masked_map(const IMP::algebra::Gaussian3Ds &gmm,
201  const Floats &weights,
202  IMP::em::DensityMap *densitymap,
203  double threshold) {
204  DensityGrid mask = IMP::bayesianem::get_grid(densitymap);
205 
206  IMP_FOREACH(const DensityGrid::Index & i, mask.get_all_indexes()) {
207  mask[i] = 0.0;
208  }
209 
210  for (unsigned int ng = 0; ng < gmm.size(); ng++) {
211  Eigen::Matrix3d covar = get_covariance(gmm[ng]);
212  Eigen::Matrix3d inverse = Eigen::Matrix3d::Zero(3, 3);
213 
214  double determinant;
215  bool invertible;
216  covar.computeInverseAndDetWithCheck(inverse, determinant, invertible);
217  IMP_INTERNAL_CHECK((!invertible || determinant < 0),
218  "\n\n\n->>>>not proper matrix!!\n\n\n");
219  double pre(get_gaussian_eval_prefactor(determinant));
220  Eigen::Vector3d evals = covar.eigenvalues().real();
221  double maxeval = sqrt(evals.maxCoeff());
222  double cutoff = 2.5 * maxeval;
223  double cutoff2 = cutoff * cutoff;
224  IMP::algebra::Vector3D c = gmm[ng].get_center();
225  IMP::algebra::Vector3D lower =
226  c - IMP::algebra::Vector3D(cutoff, cutoff, cutoff);
227  IMP::algebra::Vector3D upper =
228  c + IMP::algebra::Vector3D(cutoff, cutoff, cutoff);
229  IMP::algebra::GridIndex3D lowerindex = mask.get_nearest_index(lower);
230  IMP::algebra::GridIndex3D upperindex = mask.get_nearest_index(upper);
231  Eigen::Vector3d center(c.get_data());
232  IMP_INTERNAL_CHECK(invertible, "matrix wasn't invertible! uh oh!");
234  lowerindex, upperindex, {
235  IMP::algebra::GridIndex3D i(voxel_index[0], voxel_index[1],
236  voxel_index[2]);
237  Eigen::Vector3d r(get_vec_from_center(i, mask, center));
238  if (r.squaredNorm() < cutoff2) {
239  update_value(&mask, i, r, inverse, pre, weights[ng]);
240  }
241  })
242  }
243  Pointer<IMP::em::DensityMap> maskmap = IMP::em::create_density_map(mask);
244  Pointer<IMP::em::DensityMap> ret = IMP::em::multiply(maskmap, densitymap);
245  return ret.release();
246 }
247 
248 IMP::em::DensityMap *get_sub_map(const IMP::em::DensityMap *dm,
249  const IMP::em::DensityMap *sub_gmm,
250  const IMP::em::DensityMap *gmm) {
251  const IMP::em::DensityHeader *header = sub_gmm->get_header();
252  Pointer<IMP::em::DensityMap> m_map(IMP::em::create_density_map(dm));
253  IMP::em::emreal *data1 = sub_gmm->get_data();
254  IMP::em::emreal *data2 = gmm->get_data();
255  IMP::em::emreal *new_data = m_map->get_data();
256  for (long i = 0; i < header->get_number_of_voxels(); i++) {
257  if (data2[i] != 0.0) {
258  double const w(data1[i] / data2[i]);
259  new_data[i] *= w;
260  } else {
261  new_data[i] = 0.0;
262  }
263  }
264  return m_map.release();
265 }
266 
267 double get_overlap_fast(const IMP::algebra::Gaussian3Ds &gmm,
268  const Floats &weights, IMP::em::DensityMap *dm,
269  double factor = 2.5) {
270  DensityGrid rasterized = IMP::algebra::get_rasterized_fast(
271  gmm, weights, dm->get_header()->get_spacing(), get_bounding_box(dm),
272  factor);
273  DensityGrid orig = IMP::bayesianem::get_grid(dm);
274  /*
275  double m1(0.0);
276  double m2(0.0);
277  IMP_FOREACH(const IMP::algebra::DensityGrid::Index & i,
278  orig.get_all_indexes()) {
279  IMP::algebra::Vector3D position(orig.get_center(i));
280  IMP::algebra::DensityGrid::Index
281  j(rasterized.get_nearest_index(position));
282  m1+=rasterized[j];
283  m2+=orig[i];
284  }
285  double scale = m2/m1;
286  */
287  double score(0.0);
288  int norm(0);
289  IMP_FOREACH(const DensityGrid::Index & i,
290  orig.get_all_indexes()) {
291  IMP::algebra::Vector3D position(orig.get_center(i));
292  DensityGrid::Index j(rasterized.get_nearest_index(position));
293  // rasterized[j]*=scale;
294  double x = rasterized[j] - orig[i];
295  score += x * x;
296  ++norm;
297  }
298  return score / norm;
299 }
300 
301 IMP::algebra::Vector3Ds get_overlap_binned(const IMP::algebra::Gaussian3Ds &gmm,
302  const Floats &weights, IMP::em::DensityMap *dm,
303  double factor = 2.5, int Nbin=100) {
304 
305  Floats score(Nbin,0.);
306  Floats densities(Nbin,0.);
307  Ints counts(Nbin,0);
308 
309  double max(0.0);
310  double min(std::numeric_limits<double>::max());
311  for (long i=0 ; i<dm->get_number_of_voxels() ; ++i){
312  double const rho = dm->get_value(i);
313  if (rho>0){
314  //std::cerr << "rho " << i << " " << rho << "\n";
315  }
316  if (rho>max){
317  max=rho;
318  }
319  if ((rho>0.0) && (rho<min)){
320  min=rho;
321  }
322  }
323  double const dx( (max-min)/(Nbin) );
324  //std::cerr << "min " << min << "\nmax " << max << "\ndx " << dx << "\n";
325 
326  for(int i(0) ; i<Nbin ; ++i){
327  densities[i]=min+dx*i;
328  }
329 
330  DensityGrid rasterized = IMP::algebra::get_rasterized_fast(
331  gmm, weights, dm->get_header()->get_spacing(), get_bounding_box(dm),
332  factor);
333  IMP::em::DensityMap* rasterized_map = IMP::em::create_density_map(rasterized);
334  for (long i=0 ; i<dm->get_number_of_voxels() ; ++i){
335  double const rho = dm->get_value(i);
336  if(rho>=min){
337  size_t const index = floor((rho-min)/dx);
338  //IMP::algebra::Vector3D const position = dm->get_location_by_voxel(i);
339  double const gmm_val = rasterized_map->get_value(i);
340  double const x = gmm_val - rho;
341  score[index] += x * x;
342  ++counts[index];
343  }
344  }
345  IMP::algebra::Vector3Ds ret(Nbin);
346  for(int i(0) ; i<Nbin ; ++i){
347  if(counts[i]>0){
348  score[i]/=counts[i];
349  }
350  ret[i][0]=densities[i];
351  ret[i][1]=score[i];
352  ret[i][2]=counts[i];
353  }
354  return ret;
355 }
356 
357 IMP::algebra::Rotation3D get_rotation_matrix(const IMP::algebra::Vector3D& x,
358  const IMP::algebra::Vector3D& y) {
360  return IMP::algebra::get_rotation_from_matrix(x[0], x[1], x[2], y[0], y[1], y[2],
361  z[0], z[1], z[2]);
362 }
363 
364 float sgn(double v) {
365  return ( v < 0.0 ) ? -1.0 : ( ( v > 0.0 ) ? 1.0 : 0 );
366 }
367 
368 IMP::algebra::PrincipalComponentAnalysis NormalizePCA(const IMP::algebra::PrincipalComponentAnalysis& pca,
369  const IMP::Particles& ps)
370 {
371  // prepare vector for new pcas
372  IMP::algebra::Vector3Ds newpcs;
373  // cycle on pca
374  for (unsigned int i = 0; i < 3; ++i) {
375  // get pca
376  IMP::algebra::Vector3D x = pca.get_principal_component(i);
377  // calculate the sign of the sum of the signed inner products
378  double ips = 0.0;
379  // cycle on all particles
380  for (unsigned int j = 0 ; j < ps.size() ; ++j ) {
381  // scalar product
382  IMP::algebra::Vector3D y = IMP::core::XYZ(ps[j]).get_coordinates()-pca.get_centroid();
383  double ip = x.get_scalar_product(y);
384  // increment ips
385  ips += sgn(ip) * ip * ip;
386  }
387  // check sign of ips
388  if ( ips < 0. ) x *= -1.0;
389  // store pca
390  newpcs.push_back(x);
391  }
392  // create new pcas
393  IMP::algebra::PrincipalComponentAnalysis newpcas =
394  IMP::algebra::PrincipalComponentAnalysis(newpcs, pca.get_principal_values(), pca.get_centroid());
395 
396  return newpcas;
397 }
398 
399 
400 // MB old
401 
402 IMP::algebra::Transformation3Ds PCAalign(const IMP::algebra::PrincipalComponentAnalysis& pca1,
403  const IMP::algebra::PrincipalComponentAnalysis& pca2)
404  {
405  algebra::Transformation3Ds transforms;
406 
407  IMP::algebra::Vector3D x = pca2.get_principal_component(0);
408  IMP::algebra::Vector3D y = pca2.get_principal_component(1);
409  IMP::algebra::Rotation3D r = IMP::bayesianem::get_rotation_matrix(x, y);
411  IMP::algebra::Transformation3D(r, -(r * pca2.get_centroid())));
412  IMP::algebra::Transformation3D inverse_map_trans = map_trans.get_inverse();
413 
414  // align the principal components by enumeration 6 xy choices
415 
416  IMP::algebra::Vector3D xi = pca1.get_principal_component(0);
417  IMP::algebra::Vector3D yi = pca1.get_principal_component(1);
418  IMP::algebra::Rotation3Ds rotations(1);
419  rotations[0] = IMP::bayesianem::get_rotation_matrix(xi, yi);
420  for (unsigned int k = 0; k < rotations.size(); k++) {
421  IMP::algebra::Transformation3D points_trans =
423  rotations[k], -(rotations[k] * pca1.get_centroid()));
425  inverse_map_trans * points_trans;
426  transforms.push_back(ps2dens);
427  }
428  return transforms;
429 }
430 
431 
432 /*
433 IMP::algebra::Transformation3D PCAalign(const IMP::algebra::PrincipalComponentAnalysis& pca1,
434  const IMP::algebra::PrincipalComponentAnalysis& pca2)
435 {
436  // get the two sets of pcas
437  IMP::algebra::Vector3Ds pca1v = pca1.get_principal_components();
438  IMP::algebra::Vector3Ds pca2v = pca2.get_principal_components();
439 
440  // calculate transformations
441  for (int i = 0; i < 3; i++) {
442  pca1v[i]+=pca1.get_centroid();
443  pca2v[i]+=pca2.get_centroid();
444  }
445  pca1v.push_back(pca1.get_centroid());
446  pca2v.push_back(pca2.get_centroid());
447  algebra::Transformation3D t1 = algebra::get_transformation_aligning_first_to_second(pca1v, pca2v);
448 
449  return t1;
450 }
451 */
452 
453 
454 
455 IMPBAYESIANEM_END_NAMESPACE
456 
457 #endif /* IMPBAYESIANEM_UTILITIES_H */
Simple 3D transformation class.
const DensityHeader * get_header() const
Returns a read-only pointer to the header of the map.
Definition: DensityMap.h:243
float get_spacing() const
static const double PI
the constant pi
#define IMP_CHECK_OBJECT(obj)
Perform some basic validity checks on the object for memory debugging.
Definition: check_macros.h:277
#define IMP_FUNCTION_LOG
Beginning logging for a non-member function.
Definition: log_macros.h:293
IMP::Vector< Float > Floats
Standard way to pass a bunch of Float values.
Definition: types.h:47
emreal get_value(float x, float y, float z) const
Gets the value of the voxel located at (x,y,z)
double get_squared_distance(const VectorD< D > &v1, const VectorD< D > &v2)
Compute the squared distance between two vectors.
Definition: VectorD.h:201
#define IMP_INTERNAL_CHECK(expr, message)
An assertion to check for internal errors in IMP. An IMP::ErrorException will be thrown.
Definition: check_macros.h:139
Class for storing model, its restraints, constraints, and particles.
Definition: Model.h:72
Class for handling density maps.
Definition: DensityMap.h:93
ParticlesTemp get_selected_particles(bool with_representation=true) const
Get the selected particles.
IMP::algebra::Transformation3D get_transformation_aligning_first_to_second(const Vector3DsOrXYZs0 &source, const Vector3DsOrXYZs1 &target)
Compute the rigid transform bringing the first point set to the second.
Vector3D get_vector_product(const Vector3D &p1, const Vector3D &p2)
Return the vector product (cross product) of two vectors.
Definition: Vector3D.h:31
Ints get_index(const ParticlesTemp &particles, const Subset &subset, const Subsets &excluded)
DensityMap * multiply(const DensityMap *m1, const DensityMap *m2)
Return a density map for which voxel i contains the result of m1[i]*m2[i].
Include all non-deprecated headers in IMP.isd.
BoundingBoxD< 3 > get_bounding_box(const Cone3D &g)
Definition: Cone3D.h:64
A decorator for a particle with x,y,z coordinates.
Definition: XYZ.h:30
DensityMap * create_density_map(const DensityMap *other)
Create a copy of another map.
const algebra::Vector3D & get_coordinates() const
Convert it to a vector.
Definition: XYZ.h:109
int get_nx() const
Get the number of voxels in the x dimension.
long get_number_of_voxels() const
Get the number of map voxels.
3D rotation class.
Definition: Rotation3D.h:51
int get_ny() const
Get the number of voxels in the y dimension.
long get_voxel_by_location(float x, float y, float z) const
Calculate the voxel of a given location.
DensityGrid get_grid(IMP::em::DensityMap *in)
Return a dense grid containing the voxels of the passed density map.
Rotation3D get_rotation_from_matrix(double m00, double m01, double m02, double m10, double m11, double m12, double m20, double m21, double m22)
Generate a Rotation3D object from a rotation matrix.
IMP::Vector< Transformation3D > Transformation3Ds
Eigen::Matrix3d get_covariance(const Gaussian3D &g)
Include all non-deprecated headers in IMP.em.
VectorD< 3 > Vector3D
Definition: VectorD.h:421
IMP::algebra::Transformation3D get_transformation_aligning_first_to_second(const IMP::atom::Selection &s1, const IMP::atom::Selection &s2)
Get the transformation to align two selections.
double Float
Basic floating-point value (could be float, double...)
Definition: types.h:20
#define IMP_USAGE_CHECK(expr, message)
A runtime test for incorrect usage of a class or method.
Definition: check_macros.h:168
IMP::Vector< Int > Ints
Standard way to pass a bunch of Int values.
Definition: types.h:49
Gaussian shape.
A dense grid of values.
DenseGrid3D< double > get_rasterized_fast(const Gaussian3Ds &gmm, const Floats &weights, double cell_width, const BoundingBox3D &bb, double factor=2.5)
Rasterize the Gaussians to a grid.
int get_nz() const
Get the number of voxels in the z dimension.
#define IMP_GRID3D_FOREACH_SMALLER_EXTENDED_INDEX_RANGE(grid, center, lower_corner, upper_corner, action)
Iterate over each voxel in a subset of the grid that are less than center.
Definition: grid_utility.h:120
Select hierarchy particles identified by the biological name.
Definition: Selection.h:66