IMP logo
IMP Reference Guide  develop.e004443c3b,2024/04/25
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 (size_t 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 (size_t 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 (size_t 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 (size_t i=0; i<ps1.size(); ++i) {
86  ds1.push_back(IMP::core::XYZ(ps1[i]));
87  }
88  for (size_t i=0; i<ps2.size(); ++i) {
89  ds2.push_back(IMP::core::XYZ(ps2[i]));
90  }
91 
92  return get_rmsd_of_best_population(ds1,ds2,percentage);
93 }
94 
96  const IMP::atom::Selection &s2){
97 
100 
101  IMP::core::XYZs ds1;
102  IMP::core::XYZs ds2;
103 
104  for (size_t i=0; i<ps1.size(); ++i) {
105  ds1.push_back(IMP::core::XYZ(ps1[i]));
106  }
107  for (size_t i=0; i<ps2.size(); ++i) {
108  ds2.push_back(IMP::core::XYZ(ps2[i]));
109  }
110 
112 }
113 
114 
115 
116 double get_rmsd_of_best_population(const IMP::algebra::Vector3Ds &m1,
117  const IMP::algebra::Vector3Ds &m2,
118  const double percentage){
119  std::vector<double> sq_distances(m1.size());
120  for (size_t i=0; i<m1.size(); ++i) {
121  sq_distances[i] = get_squared_distance(m1[i], m2[i]);
122  }
123  std::sort(sq_distances.begin(), sq_distances.end());
124  double sd=0.0;
125  int N=0;
126  for (int i=0; ((double) i) < (percentage*m1.size()); ++i) {
127  sd+=sq_distances[i];
128  ++N;
129  }
130  return std::sqrt(sd/N);
131 }
132 
133 double gem_score_cc(Particles model_ps, Particles density_ps){
134  Eigen::Vector3d deriv;
135  double mm_score(0.0);
136  double md_score(0.0);
137  double dd_score(0.0);
138 
139  int nm = model_ps.size();
140  int nd = density_ps.size();
141  IMP::Model *mdl = model_ps[0]->get_model();
142 
143  for(int nd1=0 ; nd1<nd ; ++nd1){
144  for(int nd2=0 ; nd2<nd ; ++nd2){
145  dd_score += IMP::isd::score_gaussian_overlap(mdl, ParticleIndexPair(density_ps[nd1]->get_index(), density_ps[nd2]->get_index()), &deriv);
146  }
147  }
148  for(int nm1=0 ; nm1<nm ; ++nm1){
149  for(int nm2=0 ; nm2<nm ; ++nm2){
150  mm_score += IMP::isd::score_gaussian_overlap(mdl, ParticleIndexPair(model_ps[nm1]->get_index(), model_ps[nm2]->get_index()), &deriv);
151  }
152  for(int nd1=0 ; nd1<nd ; ++nd1){
153  md_score += IMP::isd::score_gaussian_overlap(mdl, ParticleIndexPair(model_ps[nm1]->get_index(), density_ps[nd1]->get_index()), &deriv);
154  }
155  }
156  double cc = md_score/std::sqrt(mm_score*dd_score);
157  return cc;
158 }
159 
160 
161 namespace {
162 double get_gaussian_eval_prefactor(double determinant) {
163  return 1.0 / pow(2 * algebra::PI, 2.0 / 3.0) / std::sqrt(determinant);
164 }
165 Eigen::Vector3d get_vec_from_center(IMP::algebra::GridIndex3D i,
166  DensityGrid const &g,
167  Eigen::Vector3d const &center) {
168  IMP::algebra::Vector3D aloc = g.get_center(i);
169  Eigen::Vector3d loc(aloc[0], aloc[1], aloc[2]);
170  Eigen::Vector3d r(loc - center);
171  return r;
172 }
173 void update_value(DensityGrid *g,
174  DensityGrid::Index i, Eigen::Vector3d r,
175  Eigen::Matrix3d inverse, double pre, Float weight) {
176  double d(r.transpose() * (inverse * r));
177  double score(pre * weight * std::exp(-0.5 * (d)));
178  (*g)[i] += score;
179 }
180 }
181 
182 DensityGrid get_grid(IMP::em::DensityMap *in) {
184  IMP_CHECK_OBJECT(in);
185  DensityGrid ret(in->get_header()->get_spacing(), get_bounding_box(in));
186  IMP_USAGE_CHECK(ret.get_number_of_voxels(0) ==
187  static_cast<unsigned int>(in->get_header()->get_nx()),
188  "X voxels don't match");
189  IMP_USAGE_CHECK(ret.get_number_of_voxels(1) ==
190  static_cast<unsigned int>(in->get_header()->get_ny()),
191  "Y voxels don't match");
192  IMP_USAGE_CHECK(ret.get_number_of_voxels(2) ==
193  static_cast<unsigned int>(in->get_header()->get_nz()),
194  "Z voxels don't match");
195  for (unsigned int i = 0; i < ret.get_number_of_voxels(0); ++i) {
196  for (unsigned int j = 0; j < ret.get_number_of_voxels(1); ++j) {
197  for (unsigned int k = 0; k < ret.get_number_of_voxels(2); ++k) {
198  DensityGrid::ExtendedIndex ei(i, j, k);
199  DensityGrid::Index gi = ret.get_index(ei);
200  long vi = in->get_voxel_by_location(i, j, k);
201  ret[gi] = in->get_value(vi);
202  }
203  }
204  }
205  return ret;
206 }
207 
208 IMP::em::DensityMap *get_masked_map(const IMP::algebra::Gaussian3Ds &gmm,
209  const Floats &weights,
210  IMP::em::DensityMap *densitymap,
211  double threshold) {
212  DensityGrid mask = IMP::bayesianem::get_grid(densitymap);
213 
214  for (const DensityGrid::Index & i : mask.get_all_indexes()) {
215  mask[i] = 0.0;
216  }
217 
218  for (unsigned int ng = 0; ng < gmm.size(); ng++) {
219  Eigen::Matrix3d covar = get_covariance(gmm[ng]);
220  Eigen::Matrix3d inverse = Eigen::Matrix3d::Zero(3, 3);
221 
222  double determinant;
223  bool invertible;
224  covar.computeInverseAndDetWithCheck(inverse, determinant, invertible);
225  IMP_INTERNAL_CHECK((!invertible || determinant < 0),
226  "\n\n\n->>>>not proper matrix!!\n\n\n");
227  double pre(get_gaussian_eval_prefactor(determinant));
228  Eigen::Vector3d evals = covar.eigenvalues().real();
229  double maxeval = sqrt(evals.maxCoeff());
230  double cutoff = 2.5 * maxeval;
231  double cutoff2 = cutoff * cutoff;
232  IMP::algebra::Vector3D c = gmm[ng].get_center();
233  IMP::algebra::Vector3D lower =
234  c - IMP::algebra::Vector3D(cutoff, cutoff, cutoff);
235  IMP::algebra::Vector3D upper =
236  c + IMP::algebra::Vector3D(cutoff, cutoff, cutoff);
237  IMP::algebra::GridIndex3D lowerindex = mask.get_nearest_index(lower);
238  IMP::algebra::GridIndex3D upperindex = mask.get_nearest_index(upper);
239  Eigen::Vector3d center(c.get_data());
240  IMP_INTERNAL_CHECK(invertible, "matrix wasn't invertible! uh oh!");
242  lowerindex, upperindex, {
243  IMP::algebra::GridIndex3D i(voxel_index[0], voxel_index[1],
244  voxel_index[2]);
245  Eigen::Vector3d r(get_vec_from_center(i, mask, center));
246  if (r.squaredNorm() < cutoff2) {
247  update_value(&mask, i, r, inverse, pre, weights[ng]);
248  }
249  })
250  }
251  Pointer<IMP::em::DensityMap> maskmap = IMP::em::create_density_map(mask);
252  Pointer<IMP::em::DensityMap> ret = IMP::em::multiply(maskmap, densitymap);
253  return ret.release();
254 }
255 
256 IMP::em::DensityMap *get_sub_map(const IMP::em::DensityMap *dm,
257  const IMP::em::DensityMap *sub_gmm,
258  const IMP::em::DensityMap *gmm) {
259  const IMP::em::DensityHeader *header = sub_gmm->get_header();
260  Pointer<IMP::em::DensityMap> m_map(IMP::em::create_density_map(dm));
261  double *data1 = sub_gmm->get_data();
262  double *data2 = gmm->get_data();
263  double *new_data = m_map->get_data();
264  for (long i = 0; i < header->get_number_of_voxels(); i++) {
265  if (data2[i] != 0.0) {
266  double const w(data1[i] / data2[i]);
267  new_data[i] *= w;
268  } else {
269  new_data[i] = 0.0;
270  }
271  }
272  return m_map.release();
273 }
274 
275 double get_overlap_fast(const IMP::algebra::Gaussian3Ds &gmm,
276  const Floats &weights, IMP::em::DensityMap *dm,
277  double factor = 2.5) {
278  DensityGrid rasterized = IMP::algebra::get_rasterized_fast(
279  gmm, weights, dm->get_header()->get_spacing(), get_bounding_box(dm),
280  factor);
281  DensityGrid orig = IMP::bayesianem::get_grid(dm);
282  /*
283  double m1(0.0);
284  double m2(0.0);
285  for (const IMP::algebra::DensityGrid::Index & i :
286  orig.get_all_indexes()) {
287  IMP::algebra::Vector3D position(orig.get_center(i));
288  IMP::algebra::DensityGrid::Index
289  j(rasterized.get_nearest_index(position));
290  m1+=rasterized[j];
291  m2+=orig[i];
292  }
293  double scale = m2/m1;
294  */
295  double score(0.0);
296  int norm(0);
297  for (const DensityGrid::Index & i : orig.get_all_indexes()) {
298  IMP::algebra::Vector3D position(orig.get_center(i));
299  DensityGrid::Index j(rasterized.get_nearest_index(position));
300  // rasterized[j]*=scale;
301  double x = rasterized[j] - orig[i];
302  score += x * x;
303  ++norm;
304  }
305  return score / norm;
306 }
307 
308 IMP::algebra::Vector3Ds get_overlap_binned(const IMP::algebra::Gaussian3Ds &gmm,
309  const Floats &weights, IMP::em::DensityMap *dm,
310  double factor = 2.5, int Nbin=100) {
311 
312  Floats score(Nbin,0.);
313  Floats densities(Nbin,0.);
314  Ints counts(Nbin,0);
315 
316  double max(0.0);
317  double min(std::numeric_limits<double>::max());
318  for (long i=0 ; i<dm->get_number_of_voxels() ; ++i){
319  double const rho = dm->get_value(i);
320  if (rho>0){
321  //std::cerr << "rho " << i << " " << rho << "\n";
322  }
323  if (rho>max){
324  max=rho;
325  }
326  if ((rho>0.0) && (rho<min)){
327  min=rho;
328  }
329  }
330  double const dx( (max-min)/(Nbin) );
331  //std::cerr << "min " << min << "\nmax " << max << "\ndx " << dx << "\n";
332 
333  for(int i(0) ; i<Nbin ; ++i){
334  densities[i]=min+dx*i;
335  }
336 
337  DensityGrid rasterized = IMP::algebra::get_rasterized_fast(
338  gmm, weights, dm->get_header()->get_spacing(), get_bounding_box(dm),
339  factor);
340  IMP::em::DensityMap* rasterized_map = IMP::em::create_density_map(rasterized);
341  for (long i=0 ; i<dm->get_number_of_voxels() ; ++i){
342  double const rho = dm->get_value(i);
343  if(rho>=min){
344  size_t const index = floor((rho-min)/dx);
345  //IMP::algebra::Vector3D const position = dm->get_location_by_voxel(i);
346  double const gmm_val = rasterized_map->get_value(i);
347  double const x = gmm_val - rho;
348  score[index] += x * x;
349  ++counts[index];
350  }
351  }
352  IMP::algebra::Vector3Ds ret(Nbin);
353  for(int i(0) ; i<Nbin ; ++i){
354  if(counts[i]>0){
355  score[i]/=counts[i];
356  }
357  ret[i][0]=densities[i];
358  ret[i][1]=score[i];
359  ret[i][2]=counts[i];
360  }
361  return ret;
362 }
363 
364 IMP::algebra::Rotation3D get_rotation_matrix(const IMP::algebra::Vector3D& x,
365  const IMP::algebra::Vector3D& y) {
367  return IMP::algebra::get_rotation_from_matrix(x[0], x[1], x[2], y[0], y[1], y[2],
368  z[0], z[1], z[2]);
369 }
370 
371 float sgn(double v) {
372  return ( v < 0.0 ) ? -1.0 : ( ( v > 0.0 ) ? 1.0 : 0 );
373 }
374 
375 IMP::algebra::PrincipalComponentAnalysis NormalizePCA(const IMP::algebra::PrincipalComponentAnalysis& pca,
376  const IMP::Particles& ps)
377 {
378  // prepare vector for new pcas
379  IMP::algebra::Vector3Ds newpcs;
380  // cycle on pca
381  for (unsigned int i = 0; i < 3; ++i) {
382  // get pca
383  IMP::algebra::Vector3D x = pca.get_principal_component(i);
384  // calculate the sign of the sum of the signed inner products
385  double ips = 0.0;
386  // cycle on all particles
387  for (unsigned int j = 0 ; j < ps.size() ; ++j ) {
388  // scalar product
389  IMP::algebra::Vector3D y = IMP::core::XYZ(ps[j]).get_coordinates()-pca.get_centroid();
390  double ip = x.get_scalar_product(y);
391  // increment ips
392  ips += sgn(ip) * ip * ip;
393  }
394  // check sign of ips
395  if ( ips < 0. ) x *= -1.0;
396  // store pca
397  newpcs.push_back(x);
398  }
399  // create new pcas
400  IMP::algebra::PrincipalComponentAnalysis newpcas =
401  IMP::algebra::PrincipalComponentAnalysis(newpcs, pca.get_principal_values(), pca.get_centroid());
402 
403  return newpcas;
404 }
405 
406 
407 // MB old
408 
409 IMP::algebra::Transformation3Ds PCAalign(const IMP::algebra::PrincipalComponentAnalysis& pca1,
410  const IMP::algebra::PrincipalComponentAnalysis& pca2)
411  {
412  algebra::Transformation3Ds transforms;
413 
414  IMP::algebra::Vector3D x = pca2.get_principal_component(0);
415  IMP::algebra::Vector3D y = pca2.get_principal_component(1);
416  IMP::algebra::Rotation3D r = IMP::bayesianem::get_rotation_matrix(x, y);
418  IMP::algebra::Transformation3D(r, -(r * pca2.get_centroid())));
419  IMP::algebra::Transformation3D inverse_map_trans = map_trans.get_inverse();
420 
421  // align the principal components by enumeration 6 xy choices
422 
423  IMP::algebra::Vector3D xi = pca1.get_principal_component(0);
424  IMP::algebra::Vector3D yi = pca1.get_principal_component(1);
425  IMP::algebra::Rotation3Ds rotations(1);
426  rotations[0] = IMP::bayesianem::get_rotation_matrix(xi, yi);
427  for (unsigned int k = 0; k < rotations.size(); k++) {
428  IMP::algebra::Transformation3D points_trans =
430  rotations[k], -(rotations[k] * pca1.get_centroid()));
432  inverse_map_trans * points_trans;
433  transforms.push_back(ps2dens);
434  }
435  return transforms;
436 }
437 
438 
439 /*
440 IMP::algebra::Transformation3D PCAalign(const IMP::algebra::PrincipalComponentAnalysis& pca1,
441  const IMP::algebra::PrincipalComponentAnalysis& pca2)
442 {
443  // get the two sets of pcas
444  IMP::algebra::Vector3Ds pca1v = pca1.get_principal_components();
445  IMP::algebra::Vector3Ds pca2v = pca2.get_principal_components();
446 
447  // calculate transformations
448  for (int i = 0; i < 3; i++) {
449  pca1v[i]+=pca1.get_centroid();
450  pca2v[i]+=pca2.get_centroid();
451  }
452  pca1v.push_back(pca1.get_centroid());
453  pca2v.push_back(pca2.get_centroid());
454  algebra::Transformation3D t1 = algebra::get_transformation_aligning_first_to_second(pca1v, pca2v);
455 
456  return t1;
457 }
458 */
459 
460 
461 
462 IMPBAYESIANEM_END_NAMESPACE
463 
464 #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:245
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:46
double get_squared_distance(const VectorD< D > &v1, const VectorD< D > &v2)
Compute the squared distance between two vectors.
Definition: VectorD.h:188
#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:86
Class for handling density maps.
Definition: DensityMap.h:95
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:71
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.
Definition: DensityMap.h:330
3D rotation class.
Definition: Rotation3D.h:52
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:408
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:19
double get_value(float x, float y, float z) const
Gets the value of the voxel located at (x,y,z)
#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:48
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:70